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SUMMARY 


Geophysical  diffraction  tomography  (GDT)  is  a  high  resolution  technique  for  quantitative 
subsurface  imaging.  The  method  is  based  on  data  derived  from  the  propagation  of  scalar 
waves  from  an  array  of  source  positions  to  an  array  of  receivers.  Images  of  spatial  variations  in 
refractive  index  are  produced  by  a  procedure  which  propagates  the  received  signal  backwards 
through  a  subsurface  cross-section.  This  is  conceptually  similar  to  the  elements  of  optical 
holography. 

An  acoustic-based  GDT  imaging  procedure  has  been  implemented  for  acoustic  waves. 
The  field  instrumentation  consists  of  a  commercially-available,  impulsive  acoustic  source,  an 
array  of  29  hydrophone/preamp  pairs,  a  custom-made  data  acquisition  system,  and  portable 
personal  computer.  Data  is  acquired  in  an  offset  vertical  seismic  profiling  configuration  with 
the  receiver  array  positioned  in  a  borehole  and  the  source  successively  fired  at  fixed  intervals 
along  a  line  on  the  ground  surface  radially  outward  from  the  borehole.  The  resulting  image 
is  of  a  two-dimensional  vertical  cross-section  spanning  the  depth  interval  of  the  receiver 
array  and  horizontal  extent  which  corresponds  to  the  length  of  the  source  line.  The  personal 
computer  is  used  to  control  the  data  acquisition  system,  execute  all  the  required  signal 
processing  steps,  and  display  the  final  image.  The  entire  procedure  can  be  implemented  in 
near  real-time  in  the  field. 

To  date,  five  field  studies  have  performed  as  part  of  the  development  and  demonstration 
of  GDT.  These  studies  have  taken  place  in  a  variety  of  geologic  settings  and  have  successfully 
imaged  naturally-occurring  and  man-made  buried  features. 
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1  INTRODUCTION 


It  is  not  necessary  to  delve  deeply  into  the  scientific  literature  to  understand  that  disposal 
of  hazardous  waste  is  a  national  problem.  Public  health  and  environmental  damages  from 
hazardous  materials  being  dumped  or  improperly  disposed  of  are  reported  daily  in  the  news. 
Cleanup  of  these  dumps  or  burial  sites  is  a  high  national  priority  and  will  remain  such  for 
decades  based  on  the  magnitude  of  the  problem.  For  a  large  number  of  these  sites,  discovering 
what  is  buried  where  poses  an  extremely  difficult  technical  challenge  (Benson  et  al.  1982). 
Remote  sensing  has  become  a  primary  means  of  investigating  the  subsurface  features  at  such 
sites.  Nonintrusive  detection  methods  offer  numerous  benefits  over  the  more  conventional 
drilling  methods.  They  are  safer  for  the  workers,  do  not  present  the  risk  of  puncturing  waste 
containers  or  liners  allowing  for  further  contamination,  and  in  many  cases  are  faster,  cheaper, 
and  more  accurate. 

Techniques  such  as  ground  penetrating  radar,  seismic  refraction,  electrical  conductivity, 
and  magnetics  have  all  been  employed  for  remote  sensing  of  buried  wastes.  The  feature 
common  to  these  and,  in  fact,  all  remote  sensing  techniques  is  that  they  employ  seismic  or 
electromagnetic  wave  energy  to  deduce  perturbations  in  subsurface  wave  propagation  prop¬ 
erties  which  infer  the  existence  and/or  location  of  the  features  of  interest.  Most  remote 
sensing  techniques  currently  used  at  hazardous  waste  sites  have  been  in  existence  many 
years  and  were  originally  developed  for  resource  exploration.  These  methods  generally  offer 
low  spatial  resolution  as  well  as  requiring  considerable  interpretive  insight.  While  research 
continues,  the  traditional  remote  sensing  techniques  appear  to  be  almost  fully  exploited. 
Ground  penetrating  radar  is  a  relatively  new  geophysical  method  which  has  achieved  consid¬ 
erable  success.  This  method  offers  greater  resolution  than  most  others;  however,  it  requires 
interpretive  insight  and  its  value  is  limited  in  certain  soil  types. 

The  technique  presented  in  this  report,  geophysical  diffraction  tomography  (GDT)  (De- 
vaney  1984;  Witten  and  Long  1986;  Witten  and  Molyneux  1988;  Witten  and  King  1988; 
Witten  and  King  1989;  King  et  al.  1989;  King  and  Witten  1989)  is  the  next  generation  of 
remote  sensing  technology.  It  begins  with  the  basic  concept  of  all  remote  sensing  in  the  use 
of  energy  waves  to  probe  the  subsurface  environment;  however,  it  progresses  beyond  existing 
methods  by  making  a  coherent  analyses  of  the  collected  signal. 

This  report  presents  an  overview  of  geophysical  tomography,  in  general,  and  a  complete 
description  of  the  physical  and  mathematical  principles  of  geophysical  diffraction  tomogra¬ 
phy.  A  complete  description  of  the  instrumentation  and  field  implementation  of  geophysical 
tomography,  along  with  a  presentation  of  field  results  to  date,  is  also  included.  The  report 
concludes  with  discussion  of  potential  applications  of  this  technique  as  weU  as  directions  for 
future  development  efforts. 
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2  BACKGROUND 

The  technique  applied  here  to  shallow,  subsurface  imaging  is  known  as  geophysical  diffrac¬ 
tion  tomography  which  is  a  generalization  of  the  reconstruction  algorithms  of  straight-ray 
tomography  (Kak  1979)  as  used  in  x-ray  CT  scanners  for  diagnostic  medicine.  Earlier  meth¬ 
ods  used  in  geophysical  tomography  were  based  on  such  straight-ray  methods  (Anderson 
and  Dziewonski  1984;  Dines  and  Lytle  1979;  Fisk  et  al.  1987).  In  order  to  better  understand 
diffraction  tomography  and  to  establish  the  motivation  for  pursuing  the  more  complex  signal 
processing  requirements  of  this  method,  it  is  enlightening  to  review  the  basis  and  limitations 
of  the  straight  ray  methods. 

The  steps  associated  with  straight- ray  image  formation  are  illustrated  in  Fig.  1. 
Reconstructions  are  developed  from  a  system  of  experiments,  each  forming  a  partial  im¬ 
age.  In  the  simple  example  presented  here,  we  will  consider  a  light  as  our  energy  source 
and  a  white  wall  as  a  detector  plane.  The  source  and  detector  maintain  a  fixed  orientation 
relative  to  each  other  but  free  to  rotate  about  the  target  which  is  taken  to  be  an  opaque  disk. 
Each  experiment  is  performed  at  a  different  instrument  orientation  relative  to  the  target.  In 
the  first  experiment,  the  upper  left  portion  of  Fig.  1,  the  detector  line  is  vertical.  The  pres¬ 
ence  of  the  target  blocks  the  light  causing  a  shadow  to  be  cast  on  the  wall.  The  first  partial 
image  is  formed  by  tracing  a  line  from  each  edge  of  the  shadow  back  to  the  source  defining 
a  triangle  as  shown  in  the  upper  right  portion  of  Fig.  1.  It  is  clear  that,  on  the  basis  of  a 
single  experiment,  the  image  is  quite  poor.  The  image  quality  can  be  significantly  improved 
by  incorporating  information  derived  from  additional  experiments.  The  lower  left  portion 
of  Fig.  1  illustrates  the  measurement  configuration  for  a  second  experiment  in  which  the 
instrumentation  has  been  rotated  90°  relative  to  the  first  experiment.  The  procedure  follows 
that  of  the  first  experiment  with  the  formation  of  a  second  partial  image;  also  a  triangle.  A 
full  reconstruction  based  on  these  two  experiments  is  formed  by  taking  the  intersection  of 
the  two  partial  images  as  shown  in  Fig.  1.  Notice  that  the  image  has  improved  significantly 
over  that  from  a  single  experiment  and  it  is  easy  to  imagine  the  increased  image  quality  that 
can  be  obtained  from  many  experiments. 

Before  considering  the  problem  of  geophysical  tomography,  it  is  worthwhile  to  point  out 
that  this  technique  is  frequently  referred  to  as  “backprojection”  since  it  projects  information 
from  the  detector  plane  back,  along  straight  lines,  to  the  source.  As  will  become  evident 
later,  this  is  different  from  diffraction  tomography  which  is  a  “backpropagation”  method. 
In  general,  backprojection  considers  semitransparent  targets  as  well  as  opaque  ones.  Rather 
than  tracing  rays  only  at  the  edges  of  the  shadow,  a  continuum  of  rays  are  traced  from  the 
entire  detector  length  back  to  the  source.  Each  ray  is  assigned  a  fixed  value  of  gray  level 
equal  to  that  measured  at  a  fixed  position  on  the  detector  plane.  For  example,  every  ray 
not  passing  through  the  target  is  assigned  a  level  of  white,  while  rays  passing  through  the 
target  each  have  an  appropriate  gray  level  which  would  only  be  black  if  the  target  is  perfectly 
opaque.  In  this  manner  CT  scanners  provide  quantitative  images  of  local  attenuation. 

In  medical  CT  scanners,  images  are  formed  based  on  measurements  of  x-ray  intensity.  For 
geophysical  measurements,  images  can  be  derived  from  straight-ray-path  algorithms  based 
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on  either  the  time  of  first  signal  arrival  or  the  amplitude  of  this  arriving  signal.  Figure  2  is 
an  example  of  a  typical  recorded  signal  as  a  function  of  time  from  a  geophysical  experiment 
with  the  time  and  amplitude  of  the  first  arrival  denoted  by  r  and  A,  respectively.  For  an 
assumed  straight  ray  path  of  length  L  as  shown  in  Fig.  3,  the  measured  amplitude  of  the 
first  arrival  A  is  related  to  the  local,  spatially-varying  attenuation  per  unit  length  a{x)  by 

A=  a{()d£  (1) 

Jo 

where  i  is  the  position  along  the  ray.  Similarly,  the  time  of  first  arrival  r  is  related  to 
local  variations  in  wave  speed  c(x)  by 

^  =  in  P) 

Local  values  of  a  and  c  can  be  reconstructed  by  an  inversion  of  eqs.  (1)  and  (2),  respectively. 
This  can  be  accomplished  by  approximating  the  continuum  of  values  of  a  and  c  by 
discrete  values  Uij  or  c,j  averaged  over  cells  defined  by  superimposing  a  grid  onto  the 
study  region.  Then  by  considering  many  ray  paths  resulting  from  a  multiplicity  of  sources 
and  receivers  and  by  approximating  the  integrals  in  eqs.  (1)  and  (2)  as  finite  sums,  a  system 
of  algebraic  equations  is  formed  which  can  be  inverted  for  each  value  of  or  c,j  by 
standard  numerical  methods.  An  alternative  inversion  scheme  can  be  applied  by  expressing 
the  integrals  in  eqs.  (1)  and  (2)  in  a  modified  form  and  then  utilizing  methods  of  Fourier 
deconvolution. 

Straight  ray  geophysical  tomography  has  successfully  been  employed  for  a  number  of 
problems  associated  with  shallow  and  deep  subsurface  investigations.  These  include  the 
mapping  of  karst  features  (Fisk  et  al.  1987),  steam  flood  zones  associated  with  oil  recovery 
(Witterholt  et  al.  1981)  and  inhomogeneities  in  the  earth’s  mantle  (Anderson  and  Dziewonski 
1984),  as  well  as  the  location  of  tunnels  and  voids  (Olehoeft  1989).  These  studies  have 
been  performed  in  an  offset  vertical  seismic  profiling  configuration  (sources  on  the  ground 
surface  and  receivers  in  a  borehole  or  vice  versa)  such  that  shown  in  Fig.  3  using  seismic 
(acoustic)  energy  sources;  and  in  a  cross- borehole  configuration  (sources  and  receivers  in 
parallel  boreholes)  as  shown  in  Fig.  4  for  both  seismic  and  radar  frequencies  electromagnetic 
sources.  Figure  5  is  an  image  of  a  tunnel  cross-section  located  in  rock  obtained  from  cross¬ 
borehole  seismic  data  and  a  backprojection  tomography  algorithm. 

While  backprojection  (straight-ray)  algorithms  of  geophysical  tomography  have,  in  some 
cases,  provided  excellent  results,  they  have  failed  in  others  due  to  limitations  in  the  method 
which  are  well  understood.  The  first  difficulty  is  that,  unlike  medical  CT  scanners  which  are 
free  to  rotate  360°  about  the  target,  geophysical  measurements  are  necessarily  constrained 
to  certain  fixed  measurement  geometries  (Figs.  3  and  4).  This  is  illustrated  in  Fig.  6  for  an 
offset  VSP  geometry.  Shown  here  is  the  approximate  image  (cross-hatched  area)  of  a  target 
(black)  obtained  from  two  experiments  (source  positions).  Note  that  the  image  here  is  quite 
poor  and  considerably  elongated  on  a  diagonal.  Also  note  that  no  substantial  improvement 
in  image  quality  could  be  obtained  by  considering  additional  source  positions  and  that  the 
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Figure  2;  Plot  of  amplitude  vs.  time  showing  amplitude  A  of  first  arriving  signal  r. 
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SOURCES  ON  THE  GROUND  SURFACE 


Figure  3:  Illustration  of  straight-ray  geophysical  tomography  concept  in  an  offset  VSP  con¬ 
figuration. 
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Figure  4:  Cross-borehole  data  acquisition  geometry. 
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Figure  5;  Image  of  tunnel  in  rock  obtained  from  the  application  of  a  straight-ray  tomography 
algorithm  in  a  cross-borehole  configuration. 


Figure  6:  Illustration  of  the  loss  of  image  resolution  associated  with  a  fixed  measurement 
configuration. 
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image  would  be  worse  if  the  target  were  deeper  and/or  further  from  the  receiver  array.  This 
difficiency  is  associated  with  the  limited  range  of  view  angles  (perspectives)  which  result 
from  a  fixed  measurement  geometry. 

A  second  and  more  fundamental  problem  associated  with  backprojection  is  that  rays  do 
not  travel  in  straight  lines  in  all  but  the  most  homogeneous  settings.  When  a  ray  encounters 
a  region  of  different  mechanical  or  electromagnetic  properties  it  will  both  bend  (refract)  on 
transmission  and  reflect  from  the  boundary  of  the  inhomogeneity.  The  degree  of  refraction 
and  the  proportion  of  incident  energy  transmitted  and  reflected  will  depend  on  the  local 
refractive  index  of  the  inhomogeneity.  The  amplitude  of  the  received  signal  will  be  reduced 
due  to  the  portion  of  the  incident  energy  which  is  reflected  and  this  will  produce  an  arti¬ 
ficially  high  value  of  imaged  local  attenuation.  Furthermore,  the  ray  bending  which  occurs 
during  transmission  will  alias  the  “shadow”  cast  by  the  target  and,  as  a  result,  it  will  be 
backprojected  along  erroneous  straight-ray  paths.  This  is  illustrated  in  Fig.  7.  Because  of 
these  limitations  of  straight  ray  algorithms,  they  are  best  suited  for  applications  in  rela¬ 
tively  simple  settings  containing  relatively  large  targets  where  high  spatial  resolution  is  not 
required  (Dines  and  Lytle  1979).  In  more  complex  settings,  some  features  of  interest  may 
be  unresolved. 
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Figure  7;  Illustration  of  ray  bending  (solid  lines)  which  occurs  as  a  result  of  changes  in 
wave-speed  and  the  resulting  erroneous  backprojection  (dashed  lines). 


3  GEOPHYSICAL  DIFFRACTION 
TOMOGRAPHY 

As  noted  earlier,  geophysical  diffraction  tomography  (GDT)  is  considered  a  backpropagation 
technique  since  it  propagates  waves  back  from  the  receivers  to  the  source  rather  than  using 
straight-ray  path  backprojection.  This  approach  is  based  on  the  scalar  wave  equation 

\^^p  +  ky{x)p  =  f{x,u:)  (3) 

where  p  is  the  pressure  (assuming  a  seismic  source)  which  has  been  Fourier  transformed  in 
time 

pU,^)  -  J  P  ,  (4) 

u!  is  the  angular  frequency,  ko  is  the  wavenumber  at  frequency  lo  and  reference  wave  speed 
Coi  kg  =  n{x)  is  the  local  value  of  refractive  index  n  =  Co(c(^),  and  /  characterizes 

the  spatial  and  frequency  distribution  of  the  wave  source.  The  objective  of  GDT  is  to  invert 
eq.  (3)  in  order  to  reconstruct  n{x)  from  values  of  the  pressure  p  measured  over  some 
receiver  contour.  This  is  accomplished  by  writing  the  formal  solution  to  eq.  (3)  (Morse  and 
Feshbach  1953)  as 

p(i,a;)  =  J  G{x-  £)  /(x,w)  dl-kl  J  Gix  -  £)  0(£)  p(£,w)  ,  (5) 

where  0{z)  =  1  —  n^{s.)  is  referred  to  as  the  “object  profile”  and  G{x)  is  the  Green’s 
function  for  propagation. 

Before  continuing  the  discussion  of  GDT  it  is  worthwhile  to  examine  eq.  (5)  and  compare 
it  to  the  backprojection  formulation  of  eqs.  (1)  and  (2).  First  consider  the  left  side  of  eq.  (5) 
which  is  composed  of  two  terms.  The  first  is  the  measured  pressure  p  and  the  second  is 
a  quantity  known  from  /  and  G.  These  terms  can  be  collectively  called  D(x, u;),  the 
reduced  data,  and  eq.  (5)  can  be  written  as 

D{x,uj)  =  -  kl  j  G{x  -  0  O(^)  p{i,  u)  di  .  (6) 

We  immediately  can  see  that  this  equation  is  quite  similar  to  those  of  backprojection,  eqs.  (1) 
and  (2),  in  that  both  formulations  have  a  known  quantity  derived  from  measurements  on  the 
left  sides  of  the  equations,  while  the  functions  of  interest  (  0(x)  for  backpropagation;  a{x) 
and  c{x)  for  backprojection)  are  under  the  integral  on  the  right  side.  Equation  (6)  is,  in 
fact,  a  generalization  of  backprojection  and  eqs.  (1)  and  (2)  can  be  derived  from  eq.  (6)  by 
invoking  the  appropriate  simplifying  assumption.  Furthermore,  for  nonattenuating  media 
0{x)  is  real;  however,  for  attenuating  media  0(x)  is  complex  with  the  real  part  related  to 
speed  variations  and  the  imaginary  part  related  to  local  variations  in  attenuation. 

One  complicating  aspect  of  backpropagation  is  that  the  pressure  appears  under  the  in¬ 
tegral  of  eq.  (6)  along  with  the  object  profile.  This  means  that  the  pressure  must  be  known 
everywhere  that  values  of  0{x)  are  desired,  not  just  along  some  measurement  contour. 
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While  it  is  possible  to  invert  this  nonlinear  formulation,  it  is  quite  difficult  as  weU  as  compu¬ 
tationally  intensive.  For  this  reason,  one  simplifying  assumption  is  made  in  the  development 
of  GDT,  the  weak  scatter  approximation.  It  is  assumed  that  the  total  pressure  p  is  com¬ 
posed  of  an  incident  pressure,  p,-,  that  satisfies  eq.  (3)  with  n{x)  =  1  (0(s.)  =  0)  along  with 
a  perturbed  field  associated  with  inhomogeneities  (nonzero  values  of  0(2.)).  It  is  further 
assumed  that  inhomogeneities  produce  a  scattered  field  which  is  small  compared  to  p,.  With 
this,  the  linearized  formulation  is 

D{x,u)  ^  -  kl  j  G{x  -  0  O(^)  p,(0  ,  (7) 

where  p,-  satisfies 

VV  +  =  /(£,^)  •  (8) 

In  this  form,  the  only  unknown  on  the  right  side  of  eq.  (7)  is  the  quantity  to  be  imaged, 

0(2). 

There  are  two  forms  of  weak  scatter  approximations  that  can  be  invoked,  the  Born  (Wolf 
1969)  and  the  Rytov  (Devaney  1981).  The  Born  approximation  assumes  that  the  total 
pressure  is  the  sum  of  the  incident  field,  pi,  and  a  scattered  field,  ps,  so  that 

P  =  P.  +  Ps  •  (9) 

The  Rytov  approximation  takes  the  perturbations  in  pressure  associated  with  inhomo¬ 
geneities  to  be  a  multiplicative  correction  to  the  incident  field  as 

p  =  p.e’^'  ,  (10) 


where  ip  is  referred  to  as  the  complex  phase.  The  imaginary  part  of  ip  represents  pertur¬ 
bations  in  the  phase  of  the  received  signal  associated  with  the  presence  of  inhomogeneities. 
While  the  real  part  of  ip  contains  similar  information  associated  with  perturbations  in 
signal  amplitude,  results  will  differ  according  to  which  weak  scatter  approximation  is  used. 
The  relative  merits  of  each  and  the  motivation  for  selecting  one  over  the  other  is  addressed 
in  the  next  chapter.  The  remainder  of  the  derivation  considered  here  can  be  cast  in  a  unified 
form  by  the  definition  of  the  reduced  data,  D,  given  by 


D{x,io)  = 


Ps  ,  Born  approximation] 
Piip  ,  Rytov  approximation. 


(11) 


The  inversion  of  eq.  (7)  is  straight  forward  but  tedious.  Therefore,  the  detailed  mathe¬ 
matics  for  relatively  general  measurement  geometries  is  provided  in  Appendix  A.  Here  we 
highlight  the  inversion  procedure  for  an  offset  VSP  configuration  focusing  on  the  physical 
interpretation  of  the  various  steps  in  the  procedure. 

Conceptually,  geophysical  diffraction  tomography  is  best  described  as  an  analog  to  optical 
holography  as  depicted  in  Fig.  8  (King  et  al.  1989).  The  first  step  in  the  imaging  procedure 
is  to  construct  the  reduced  data  D  [eq.  (11)]-  This  is  analogous  to  the  interference  pattern 
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Figure  8:  Illustration  of  the 


physical  concept  of  geophysical  diffraction  tomography. 
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RECEIVER  ARRAY 


created  in  optical  holography  by  the  interaction  between  the  laser  illuminating  the  target 
and  a  reference  beam.  The  next  two  mathematical  steps  associated  with  the  inversion  of 
eq.  (7)  are 

7  (<'  .  io)  =  y  4)  ■  a/  D(e  ,  ,  (12) 

the  synthetic  aperture  step  (Schultz  and  Clarebout  1978);  and 

d[K{s  -  ^)]  =  ^  m  •  6  /  dt  ,  (13) 

Kq  J 

the  holographic  lens  step.  In  eq.  (12)  and  (13),  to  and  f  are  the  positions  of  sources  and 
receivers,  respectively;  and  other  parameters  are  defined  in  Fig.  9.  The  step  corresponding 
to  eq.  (12)  is  a  coherent  sum  or  synthetic  aperture  step  synthesizing  a  coherent  wave  from  a 
superposition  of  all  sources.  Thus,  this  step  produces  an  incident  wave  form  similar  to  the 
coherent  light  produced  by  laser  illumination.  Different  viewing  perspectives  are  obtained  by 
varying  the  direction  of  propagation  ^  of  this  synthetic,  illuminating  incident  wave  field. 
Equation  (13)  is  the  mathematical  analog  of  a  holographic  lens,  refocusing  the  wavefield 
distorted  by  inhomogeneities.  The  result  of  the  application  of  eq.  (13)  is  the  spatial  Fourier 
transform  0{K)  of  the  object  profile  0{x).  This  is  similar  to  the  image  contained  in  the 
optical  plate  of  holography.  Finally,  the  image  can  be  recovered  by  numerical  inversion. 

0{x)  =  j  d{K)  •  ^dK  ,  (14) 

which  is  the  mathematical  analog  to  the  presentation  of  the  hologram  by  reillumination  of 
the  optical  plate. 

At  this  point  it  is  appropriate  to  evaluate  GDT  in  light  of  the  problems  and  limitations 
associated  with  straight-ray  backprojection  described  in  the  previous  chapter.  First  and 
foremost  are  the  complications  arising  from  ray  bending.  Since  GDT  is  a  backpropagation 
method,  in  contrast  to  backprojection,  it  rigorously  accounts  for  ray  bending  within  the 
context  of  the  invoked  weak  scatter  approximation.  Thus,  no  image  artifacts  occur  provided 
that,  in  principle,  the  conditions  of  the  weak  scatter  approximation  are  satisfied.  In  practice, 
this  means  that  ray  bending  phenomenon  are  properly  treated  provided  that  strong  scatters 
do  not  occur  in  close  proximity  to  one  another.  In  situations  where  this  is  not  the  case, 
proximate  strong  scatterers  (inhomogenities)  can  be  expected  to  blur  together  and  a  weak 
scatterer  may  not  properly  be  resolved  if  it  is  proximate  to  a  strong  scatterer.  While  the 
limited  view  angles  still  occur  in  backpropagation,  the  inherent  focus  step,  eq.  (13),  provides 
better  resolution  per  experiment  than  does  backprojection  (Devaney  1983).  Rather  than  a 
triangular  partial  reconstruction  as  depicted  in  Figs.  1  and  6,  a  single  partial  image  of  a 
circular  disk  derived  from  GDT  will  be  eliptic. 

One  final  and  important  point  must  be  mentioned  when  comparing  GDT  with  straight- 
ray  tomography.  Unlike  the  straight-ray  approach  which  is  based  on  the  first  arriving  signal, 
GDT  is  a  full  waveform  technique  extracting  information  from  the  entire  time  series  (Fig.  2). 
It  is  clear  from  eq.  (4)  that  all  features  in  the  time  series  which  influence  the  pressure  at 
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Figure  9:  Notation  used  in  the  formulation  of  geophysical  diffraction  tomography  in  an  offset 
VSP  configuration. 
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SOURCES  ON  THE  GROUND  SURFACE 


Figure  10;  Illustration  of  possible  multiple  ray  paths  which  produce  constructive  and/or 
destructive  interference  in  the  recorded  time  series. 


RECEIVERS  IN  A  BOREHOLE 


a  fixed  frequency  uj  are  considered  in  the  GDT  algorithm.  This  is  illustrated  in  Fig.  10. 
Shown  here  are  two  realizable  ray  paths;  one  a  direct  path  and  the  other  refracted  through  the 
inhomogeneity.  Under  most  conditions,  the  direct  ray  will  be  that  which  produces  the  first 
arriving  signal.  Consequently,  this  path  provides  little  information  about  the  inhomogeneity 
(except  where  it  isn’t)  when  considering  first  arrival  information.  In  contrast,  GDT  extracts 
information  for  both  ray  paths  as  well  as  many  others.  Furthermore,  unlike  straight-ray 
tomography,  GDT  is  not  strictly  a  transmission  technique.  Therefore  supplemental  view 
angles  can  be  considered  by  inclusion  of  slant-stack  or  synthetic  aperture  directions  which 
are  away  from  the  receiver  array.  As  illustrated  in  Fig.  11,  this  allows  the  incorporation  of 
other  late  arriving  signals  associated  with  reflections  or  back  scatter. 

The  signal  processing  steps  associated  with  the  implementation  of  the  GDT  algorithm 
described  above  are  summarized  in  Fig.  12. 
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Figure  12:  Sequence  of  steps  employed  in  the  implementation  of  geophysical  diffraction 
tomography. 
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4  IMPLEMENTATION  SCHEMES, 
INSTRUMENTATION,  AND  FIELD 
OPERATIONS 

4.1  Measurement  Geometries 

There  are  three  basic  practical  instrument  configuration  for  field  implementation  of  geophys¬ 
ical  diffraction  tomography.  These  are  cross-borehole,  offset  VSP,  and  surface-to-surface  ge¬ 
ometries.  The  surface-to-surface  configuration  employs  source  and  receiver  arrays  which  are 
both  located  on  the  ground  surface  [Fig.  13(a)].  In  the  cross-borehole  configuration,  an  array 
of  acoustic  sources  is  located  in  one  borehole  and  an  array  of  acoustic  receivers  is  placed  in 
a  parallel  borehole  [Fig.  13(b)].  For  offset  VSP  [Fig.  13(c)],  receivers  are  again  located  in  a 
borehole  but  sources  are  deployed  on  a  line  on  the  ground  surface  extending  outward  from 
the  borehole. 

The  ultimate  regions  which  can  be  imaged  are  defined  by  the  particular  geometry  and 
the  lengths  of  both  the  source  and  receiver  arrays.  For  the  cross-borehole  geometry  this  two- 
dimensional  vertical  cross-section  is  a  rectangular  region  with  a  vertical  extent  corresponding 
to  the  depth  interval  spanned  by  the  receiver  array  and  a  horizontal  extent  defined  by  the 
distance  between  boreholes.  For  the  offset  VSP  configuration,  the  vertical  extent  of  the  image 
is  the  depth  interval  spanned  by  the  receivers  while  the  horizontal  boundaries  of  the  imaged 
region  are  defined  by  the  length  of  the  source  array  on  the  ground  surface.  The  vertical 
extent  of  the  region  which  can  be  imaged  in  the  surface-to-surface  geometry  is  arbitrary  and 
the  horizontal  extent  of  this  region  corresponds  to  the  interval  spanned  by  the  surface-located 
receivers.  In  practice,  however,  these  regions  could  be  somewhat  smaller  due  to  the  finite 
length  of  the  source  array.  This  is  because  it  is  necessary  to  simulate  plane  wave  illumination 
by  means  of  a  slant  stack  or  synthetic  aperture  procedure.  As  illustrated  in  Fig.  14,  this 
synthetic  wave  is  only  near-planar  over  a  certain  portion  of  the  illuminated  region.  It  is 
therefore  necessary  to  either  reduce  the  imaged  region  or  restrict  the  range  of  viewing  angles 
in  order  to  minimize  image  artifacts. 

While  each  configuration  described  above  can  ideally  achieve  a  range  of  viewing  angles 
approaching  180°,  there  exists  significant  differences  in  image  quality  among  these  mea.sure- 
ment  geometry.  This  is  illustrated  in  Fig.  15  which  shows  images  of  a  single  circular  disk  for 
a  number  of  geometries.  The  reason  for  these  evident  differences  is  that,  while  each  config¬ 
uration  can  utilize  a  comparable  range  of  view  angles,  the  specific  view  angles  vary  among 
configurations.  It  is  clear  that  the  best  image  is  obtained  from  a  cross- borehole  configura¬ 
tion  [Fig.  15(a)]  while  the  worst  is  offered  by  the  surface-to-surface  method  [Fig.  15(c)].  This 
difference  is  because  the  quality  of  information  contained  in  a  single  view  angle  is  a  function 
of  the  particular  angle.  The  most  useful  information  is  contained  in  transmission  angles 
while  the  value  of  reflection  or  backscatter  angles  is  minimal.  For  cross-hole  geometries  all 
available  view  angles  are  associated  with  transmission  and,  for  offset  VSP,  50%  of  available 
views  are  transmission  angles.  The  surface-to-surface  configuration  offers  no  transmission 
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Figure  13:  Notation  used  in  formulating  the  geophysical  diffraction  tomography  algorithm 
in  (a)  surface-to-surface,  (b)  cross- borehole,  and  (c)  offset  VSP  geometry. 
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Figure  14;  “Edge”  effects  cissociated  with  the  synthesis  of  a  plane  wave  by  slant-stacking. 
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Figure  15:  Synthetic  images  of  a  circular  disk  for  a  (a)  cross-borehole,  (b)  offset  VSP. 
(c)  surface-to-surface,  and  (d)  composite  cross- borehole  and  offset  VSP  geometry. 
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angles.  It  is  clear  that  the  differences  in  image  quality  between  cross- borehole  [Fig.  15(a)] 
and  offset  VSP  [Fig.  15(b)]  are  slight;  however,  the  image  quality  attainable  from  purely  sur¬ 
face  measurements  [Fig.  15(c)]  is  unsuitable  for  most  applications.  In  summary,  reflection 
view  angles  are  of  little  value  as  a  primary  source  of  information  and  are  best  utilized  as  a 
supplement  to  sharpen  images. 

Along  with  the  three  primary  configurations  described  above,  composite  geometries  can 
be  used  to  further  improve  image  quality.  One  such  composite  geometry  is  surface  to  two 
borehole  which  provides  the  same  range  of  view  angles  as  offset  VSP;  however,  here  improved 
image  quality  results  because  all  view  angles  are  associated  with  transmission.  An  example 
of  the  use  of  this  geometry  is  given  in  the  next  chapter.  Another  composite  system  is 
combined  cross-hole  and  offset  VSP  where  180°  of  transmission  angles  and  an  additional  90° 
of  backscatter  angles  can  be  realized.  Figure  15(d)  is  an  example  of  a  simulated  image  from 
this  geometry. 

The  offset  VSP  configuration  appears  to  be  best  suited  for  shallow  applications  associated 
with  environmental  problems  because  it  provides  good  image  quality,  avoids  the  complica¬ 
tions  of  a  downhole  energy  source,  and  allows  for  the  imaging  of  multiple  cross-sections  from 
a  single  borehole.  For  imaging  in  a  cross-borehole  configuration,  two  boreholes  are  required 
to  image  each  cross-section  and  one  additional  borehole  would  be  required  for  each  additional 
cross-section.  Multiple  cross-sections  or  a  three-dimensional  image  can  be  obtained  in  the 
offset  VSP  configuration  by  the  deployment  scheme  illustrated  in  Fig.  16. 

The  remainder  of  this  chapter  deals  with  the  field  instrumentation  and  implementation 
for  GDT  in  an  offset  VSP  configuration. 

4.2  Instrumentation 

The  data  acquisition  system  consists  of  the  four  hardware  subsystems  represented  in  Fig¬ 
ure  17  plus  the  associated  operating  software.  Each  of  the  separate  components  are  described 
below. 

Noise  Source  —  The  propagating  sound  wave  is  produced  by  a  commercially  manufac¬ 
tured  8-gauge  seismic  gun  known  as  “BETSY.”  BETSY,  shown  in  Fig.  18  fires  a  85  g  slug 
downward  producing  a  hemi-spherical  wave  front  with  an  initial  energy  of  12.2  KJoule.  The 
normal  frequency  spectrum  for  this  source  shows  the  majority  of  the  energy  is  in  the  25-250 
Hz  range.  A  typical  signature  from  one  source  firing  is  shown  in  Fig.  19.  Two  critical  param¬ 
eters  associated  with  seismic  gun  operations  are  proper  coupling  to  the  ground  and  accurate 
location  of  the  source  with  respect  to  the  borehole.  Coupling  determines  the  efficiency  in 
transferring  energy  to  the  ground.  This  directly  impacts  the  signal-to-noise  ratio  and  the 
effective  distances  for  wave  propagation.  Location  of  the  source  is  an  important  parame¬ 
ter  affecting  variability  of  results.  Sampling  procedures  present  in  the  following  subsection 
describe  methods  utilized  to  properly  locate  and  couple  the  source. 

Hydrophone  Receivers  —  Signal  collection  is  accomplished  by  29  hydrophones  and  pream¬ 
plifiers  assembled  into  a  sealed,  oil-filled  streamer  (Fig.  20).  These  hydrophones  are  com¬ 
mercially  available  and  offer  good  frequency  response  in  the  range  from  several  Hz  to  several 
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17:  Schematic  of  the  field  system  used  to  implement  geophysical  diffraction  tomog- 
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Figure  19:  Typical  time  series  and  power  spectrum  produced  by  the  BETSY  seisgun  (Source 
Martin,  P.  N.,  “Betsy  M3  seisgun  source,”  Betsy  Seisgun  Inc.,  October  1986). 
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Figure  20;  Close-up  photograph  of  a  portion  of  the  hydrophone/preamp  streamer. 
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KHz.  Each  set  is  hardwired  to  a  separate  channel  of  the  data  acquisition  system.  The 
streamer  is  lowered  into  a  water-filled  hole  to  the  desired  depth  (Fig.  21).  Utilizing  a  water- 
filled  hole  provides  a  stronger  received  signal  by  promoting  better  acoustic  coupling  between 
hydrophones  and  the  formation. 

Data  Acquisition  System  DAS  —  The  DAS  consists  of  32  identical  data  collection  chan¬ 
nels  assembled  on  8  separate  computer  boards.  Each  channel  includes  an  analog-to-digital 
converter,  an  amplifer,  local  memory,  and  access  to  a  single  board  computer.  Channels  are 
programmable  for  signal  gain,  sampling  period,  and  number  of  sample  points.  Data  collec¬ 
tion  is  initiated  by  the  system’s  electronic  trigger  actuated  by  a  signal  from  an  accelerometer 
attached  to  the  BETSY.  Figure  22  shows  the  DAS  along  with  the  supervisory  computer 
being  operated  in  the  field. 

Supervisory  Computer  —  Because  the  hardware  system  was  necessarily  constructed  be¬ 
fore  final  field  methods  were  known,  the  system  was  designed  to  provide  maximum  flexibility 
in  operation.  Key  to  obtaining  this  flexibility  was  incorporation  of  on-line  control  through  a 
supervisory  computer.  In  general,  the  computer  is  used  to  input  the  data  collection  parame¬ 
ters  executed  by  the  DAS  and  to  provide  permanent  memory  for  the  collected  data.  Current 
system  configuration  employs  a  COMPAQ  Portable  386  personal  computer  equipped  with 
a  data  acquisition  board  specially  designed  to  communicate  with  the  DAS.  The  computer 
includes  a  40  mega-byte  permanent  memory  and  one  high  density  diskette  drive.  Specifics  of 
operation  are  provided  in  the  following  section  which  described  the  DAS  controller  functions 
operating  on  the  supervisory  computer.  For  efficiency  all  data  acquisition  software  is  written 
in  FORTH,  a  relatively  low-level  programming  language. 

4.3  Field  Operations 

The  field  operations  necessary  for  implementation  of  GDT  include  site  preparation,  defi¬ 
nition  of  source  lines,  data  acquisition,  and  signal  processing.  The  elements  are  described 
individually  below. 

Site  preparation  —  It  is  necessary  to  suitably  prepare  a  site  for  the  deployment  of  sources 
and  receivers.  This  requires  a  4  in  diameter  borehole  to  accommodate  the  receiver  array.  A 
dry  or  uncased  hole  can  be  used  but,  as  noted  above,  a  water-filled  hole  is  more  desireable. 
When  preparing  a  site  from  scratch,  a  6  in  diameter  hole  is  drilled  and  a  4  inch  id  PVC 
casing  capped  on  the  bottom  is  placed  in  the  hole.  The  annular  region  between  the  casing 
on  the  borehole  is  backfilled  with  sand  or  available  site  materials  to  improve  acoustic  wave 
coupling.  It  is  not  critical  to  the  method  to  rigorously  follow  these  borehole  development 
procedures.  Suitable  monitoring  wells  can  be  used  as  well  as  open  core  holes  if  available. 

Definition  of  source  lines  —  The  acoustic  source  (BETSY)  must  be  fired  at  predetermined 
intervals  along  a  line  radially  outward  from  the  borehole  containing  the  hydrophone  array. 
While  this  line  does  not  have  to  be  horizontal  it  should  be  as  flat  as  possible.  Energy 
produced  by  the  seismic  gun  and  propagated  for  some  subsurface  distance  without  substantial 
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Figure  21:  Photography  of  the  streamer  being  lowered  into  a  monitoring  well. 
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Figure  22:  Photograph  of 


:  data  acquisition  system  (DAS)  and  supervisory  computer. 
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attenuation  losses  lies  in  the  10-300  Hz  spectral  range.  The  resolution  limits  of  GDT  are 
wavelength-dependent  with  the  smallest  inclusion  which  can  reliably  be  imaged  being  about 
a  quarter  of  a  wavelength.  The  actual  wavelengths  which  can  be  realized  are  defined  by 
the  relationship  A  =  Co//,  where  A  is  the  wavelength,  cq  is  the  reference  sound  speed 
at  the  site,  and  /  is  the  frequency  considered.  Thus,  the  best  possible  resolution  which 
can  be  achieved  is  defined  by  the  maximum  possible  frequency  which  can  be  propagated.  In 
practice,  however,  there  is  a  trade-off  because  higher  frequencies  are  more  rapidly  attenuated 
and,  as  a  result,  this  limits  the  horizontal  extent  of  the  region  which  can  be  imaged  from 
a  single  borehole.  Furthermore,  the  implementation  of  GDT  requires  the  spatial  resolution 
of  the  propagated  wave  which,  in  turn,  requires  source  positions  to  be  established  at  half 
wavelength  intervals.  In  other  words,  increased  resolution  requires  more  data  collection 
and  limits  the  distance  the  source  can  be  moved  from  the  borehole.  These  points  must  be 
considered  when  establishing  the  source  deployment  scheme. 

The  above  discussion  also  influences  site  preparation.  A  source  line  must  be  established 
with  source  firing  locations  occurring  every  half  wavelength.  Minor  deviations  in  source 
position  are  acceptable  provided  such  deviations  are  not  more  than  about  one  eighth  of  a 
wavelength.  Surface  irregularities  and  obstructions  such  as  trees  can  be  avoided  by  relocation 
of  the  source  provided  that  such  a  deviation  is  much  smaller  than  one  wavelength.  As  a  result, 
at  some  sites  minor  (hand)  grading  may  be  required  or  imaging  must  be  accomplished  at  a 
lower  frequency. 

Figure  23  is  a  photograph  of  a  prepared  site  with  cased  boreholes  in  place  and  grading 
stakes  located  to  guide  source  positioning. 

Data  acquisition  —  The  data  acquisition  system  developed  for  the  initial  field  implementation 
of  GDT  was  designed  to  be  quite  flexible  and  a  recent  modification  to  the  software  which 
controls  the  acquisition  of  data  makes  use  of  the  system  both  efficient  and  user-friendly. 
Figure  24  is  a  photograph  of  the  user  menu  appearing  on  the  monitor  of  the  supervisory 
computer.  The  shaded  area  occupying  most  of  the  screen  is  used  to  input  and  display  selected 
data  acquisition  parameter,  while  the  box  in  the  lower  right  corner  displays  supplemental 
commands.  The  elements  of  the  main  menu  and  supplement  commands  are  described  below. 

Main  Menu 

1.  File  name:  The  user  inputs  a  master  file  name  for  each  cross-section  to  be  surveyed. 
For  example  a  file  name  TEST  may  be  used.  A  separate  file  of  time  series  from  all 
hydrophones  is  created  and  stored  on  the  hard  disk.  The  procedure  assumes  that 
the  first  location  for  the  source  is  the  one  nearest  to  the  borehole.  Following  each 
data  acquisition,  data  is  stored  in  a  file  TEST01.DAT  for  the  first  source  position,  file 
TEST02.DAT  for  the  second  source  position,  etc.  Once  data  acquisition  begins,  the 
file  name  corresponding  to  a  particular  source  location  is  automatically  displayed. 

2.  No.  of  samples:  The  user  may  input  the  total  number  of  samples  in  time  (up  to  512) 
to  be  digitized  and  stored  for  each  hydrophone. 
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Figure  23:  Photograph  of  the  Chestnut  Ridge  site  as  prepared  for  the  tomography  field 
study. 
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Figure  24:  Illustration  of  the  interactive  user  menu  which  controls  data  acquisition. 
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3.  Sampling  interval:  The  user  inputs  the  time  interval  in  microseconds  between  each 
digitization.  The  total  length  of  the  acquired  time  series  (in  microseconds)  is  the 
product  of  the  sampling  interval  and  number  samples. 

4.  Initial  source  position:  The  user  inputs  the  distance  of  the  source  from  the  borehole 
at  its  initial  location. 

5.  Source  spacing:  The  user  inputs  the  distance  between  successive  source  positions. 

6.  Source  index:  Displayed  here  is  an  index  which  defines  the  incremental  source  position. 
This  is  1  for  the  first  source  position,  2  for  second  position,  etc.  The  user  may  change 
this  value  to  access  data  collected  for  a  previous  source  location.  Changing  the  source 
index  will  result  in  a  change  in  the  displayed  file  name  to  correspond  to  the  data  file 
being  accessed. 

7.  Source  position:  Displayed  here  is  the  distance  between  the  source  and  the  borehole 
for  the  current  source  position.  This  is  automatically  updated  based  on  the  current 
source  index,  initial  source  position,  and  source  spacing.  Changing  the  source  index 
or  file  name  will  result  in  a  corresponding  change  in  the  displayed  file  name,  source 
index,  and  source  position. 

8.  Frequency  point:  The  user  inputs  a  frequency  index  at  which  imaging  is  to  be  per¬ 
formed.  This  is  stored  for  use  in  later  signal  processing. 

9.  Sound  speed:  The  user  inputs  a  reference  sound  speed  for  the  site  being  studied.  This 
is  stored  for  use  in  later  signal  processing. 

10.  First  depth:  Depth  of  shallowest  hydrophone  in  the  array. 

11.  PGA  file:  Name  of  file  containing  gain  settings  for  the  programmable  gain  amplifiers. 

Supplemental  Commands 

1.  Ready:  This  command  downloads  sampling  parameters  from  the  supervisory  computer 
to  the  DAS  prior  to  source  firing.  Following  data  acquisition,  the  user  is  prompted  to 
strike  any  key.  The  response  to  this  prompt  causes  the  data  retained  in  local  memory 
of  the  DAS  to  be  downloaded  into  the  current  file  name  displayed  on  the  main  menu 
and  stored  on  the  hard  disk.  Added  to  each  data  file  is  a  heading  containing  all  unique 
parameters  currently  displayed  on  the  main  menu.  This  also  automatically  increments 
the  file  name,  source  index,  and  source  position  displayed  on  the  main  menu. 

2.  Display:  The  display  command  allows  the  user  to  view  a  plot  of  the  time  series  for 
any  hydrophone  from  the  file  currently  displayed  under  file  name  (Fig.  25).  Display  is 
typically  used  periodically  to  verify  the  system  is  acquiring  data  correctly  and  to  check 
that  the  gain  setting  is  appropriate.  The  user  may  review  data  acquired  for  a  previous 
position  by  changing  the  displayed  file  name  or  source  index. 
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3.  Unwrap:  This  command  executes  an  initial  signal  processing  step.  This  is  the  com¬ 
putation  of  the  principal  value  of  the  acoustic  wave  phase  and  its  total  value  obtained 
by  a  phase  unwrapping  procedure  (see  Appendix  A).  A  plot  of  the  unwrapped  phase 
as  a  function  of  receiver  position  appears  on  the  PC  monitor  following  completion  of 
these  computations  (Fig.  26).  The  unwrapped  phase  data  is  stored  separately  in  a  two- 
dimensional  array.  Following  each  source  firing  and  phase  unwrap,  another  column  is 
added  to  this  array.  This  processed  data  is  used  by  the  imaging  software  for  all  ac¬ 
quired  data.  The  computational  speed  of  the  COMPAQ  Portable  386  allows  the  phase 
unwrapping  to  be  accomplished  during  the  time  required  to  reposition  the  source.  Per¬ 
forming  this  operation  during  the  data  acquisition  phase  saves  about  15  min  in  signal 
processing  time  later. 

4.  Detect:  Relative  changes  in  phase  as  a  function  of  depth  between  adjacent  or  nearby 
source  positions  suggests  the  presence  of  inhomogeneities.  The  approximate  location 
and  nature  (layers,  isolated  inclusions,  etc.)  of  inhomogeneities  may  be  inferred  by 
visual  examination  of  selected  phase  plots  such  as  the  one  shown  in  Fig.  26.  By 
invoking  Detect,  the  user  can  view,  simultaneously,  up  to  six  selected  phase  plots  in 
order  to  establish  preliminary  information  regarding  subsurface  conditions.  Detect  can 
be  invoked  any  time  during  the  data  acquisition  operation.  Figure  27  is  an  example  of 
the  display  created  by  Detect. 

5.  Info:  Following  completion  of  data  acquisition  and  signal  processing  all  raw  and  pro¬ 
cessed  data  can  be  archived  on  diskette.  The  Info  command  allows  the  user  to  inter¬ 
actively  create  a  text  file  describing  experimental  conditions  and  other  observations. 
This  file  can  then  be  incorporated  in  the  information  contained  on  diskette. 

6.  Parameter  set:  This  command  calls  another  menu  which  allows  the  user  to  establish 
gain  files  for  the  programmable  gain  amplifiers  and  to  invoke  diagnostic  software  used 
to  verify  proper  communications  between  the  supervisory  PC  and  the  DAS. 

The  above  described  data  acquisition  procedure  was  designed  to  be  both  efficient  and 
user-friendly.  Data  acquisition  proceeds  with  a  minimum  of  user  intervention  so  there  is 
little  opportunity  for  user-induced  error.  The  user  has  displayed  before  him  at  all  times  all 
relevant  parameters  associated  with  the  experimental  procedure  and  appropriate  diagnostic 
functions  to  allow  easy  and  rapid  identification  of  system  malfunctions.  This  data  acquisition 
procedure  is  quite  fast  with  a  proven  speed  (including  data  acquisition,  data  storage,  and 
preliminary  signal  processing)  of  one  source  firing  at  least  every  minute. 

Signal  processing  —  Following  the  completion  of  data  acquisition  for  an  individual  surveyed 
cross-section,  the  sequence  of  signal  processing  steps  outlined  in  Fig.  28  can  be  executed  to 
create  images  in  near  real  time  in  the  field.  The  mathematical  details  associated  with  these 
steps  are  given  in  Appendix  A. 

The  first  step  of  the  signal  processing  procedure  is  to  compute  the  perturbed  phase.  This 
is  the  difference  between  the  actual  total  as  determined  by  the  phase  unwrapping  procedure 
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Figure  26:  Unwrapped  phase  as  a  function  of  receiver  index  produced  by  execution  of  the 
supplemental  UNWRAP  command. 
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Figure  27:  Example  of  graphical  output  produced  by  execution  of  the  supplemental  DETECT 
command. 
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Figure  28:  Schematic  of  the  sequence  of  signal  processing  steps  necessary  to  produce  an 
image  following  data  collection. 
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and  that  which  would  have  been  measured  if  propagation  occured  in  a  homogeneous  medium. 
Recalling  the  discussion  of  the  physical  analogy  of  GDT  to  optical  holography,  given  in  the 
previous  chapter,  this  step  is  the  one  which  corresponds  to  the  formation  of  an  interference 
pattern  created  by  the  interaction  action  of  the  target  laser  beam  with  the  reference  beam. 

The  next  step  in  the  procedure  is  the  slant-stack  to  simulate  plane  (coherent)  wave 
illumination  [eq.  (12)]  and  the  computation  of  scattering  amplitude  [the  application  of  the 
numerical  holographic  lens,  eq.  [(13)].  In  this  step  the  user  must  specify  the  number  of 
viewing  angles  to  be  used  as  well  as  specify  a  range  of  view  angles  between  0°  and  180°.  This 
selection  must  be  made  judiciously  to  avoid  image  artifacts. 

The  third  step  in  the  procedure,  is  the  numerical  inversion  of  the  spatial  Fourier  transform 
of  the  object  profile  [eq.  (14)].  The  user  specifies  the  horizontal  extent  and  resolution  of  this 
computation.  The  selection  of  the  horizontal  extent  must  be  compatible  with  the  range  of 
view  angles  selected  in  the  previous  step  to  avoid  image  artifacts. 

The  final  step  in  this  procedure  is  the  graphical  display  of  the  resulting  image.  The  means 
of  display  utilized  in  the  field  is  a  pseudo  gray  scale  image  produced  on  the  PC  monitor  with 
dark  levels  of  gray  corresponding  to  increasing  relative  values  of  sound  speed.  A  hard  copy 
of  the  result  can  be  obtained  in  the  field  through  the  use  of  a  dot  matrix  printer  and  the 
PRINT  SCREEN  command.  Figure  29  is  an  example  of  an  image  as  produced  in  the  field. 
The  image  is  somewhat  “blocky”  and  limited  to  only  15  gray  levels.  This  is  a  result  of  the 
limited  number  of  addressable  pixels  available  on  the  monitor  of  the  supervisory  PC.  The 
image  quality  can  be  improved  by  using  recently-available  portable  PC’s  with  high-resolution 
monitors.  As  will  be  seen  in  the  next  chapter,  final  images  of  far  superior  quality  can  be 
achieved  using  a  dot-matrix  or  laser  printer. 

The  signal  processing  software  described  here  is  quite  efficient,  typically  requiring  about 
2  min  for  all  data  analysis  and  display. 
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Figure  29:  Example  of  an  image  as  it  appears  on  the  monitor  of  the  supervisory  computer. 
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5  RESULTS 


This  chapter  presents  results  to  date  on  the  implementation  of  geophysical  diffraction  to¬ 
mography  from  five  field  tests  performed  at  four  sites.  These  are  presented  in  chronological 
order  and  each  successive  test  was  executed  to  either  validate  or  improve  the  technology.  The 
images  present  here,  unless  otherwise  noted,  display  gray  scales  where  increasingly  darker 
shades  of  gray  indicate  increasing  sound  speeds  of  subsurface  features.  In  interpreting  these 
images  it  is  important  to  realize  that  sound  speed  is  a  composite  function  of  both  density 
and  compressibility.  Thus,  an  elevated  sound  speed  relative  to  background  can  be  the  result 
of  an  increased  density,  decreased  compressibility,  or  both.  The  opposite  is  true  for  features 
displaying  sound  speeds  which  are  less  than  background. 

5.1  Chestnut  Ridge  March  1987  Survey 

This  field  study  was  conducted  on  the  Oak  Ridge  reservation  at  a  site  on  Chestnut  Ridge 
between  the  Y-12  and  X-10  areas.  The  target  was  a  0.6  m  diameter  cast  iron  water  pipe 
buried  approximately  1.2  m  deep  in  clay  soil.  Two  0.15  m  by  10  m  deep  boreholes  were 
installed,  one  on  each  side  of  the  cleared  pipeline  right-of-way.  PVC  casings  equipped  with 
bottom  caps  were  installed  in  each  hole  and  the  outside  of  the  wells  were  grounted  with  sand. 
The  distance  between  the  two  boreholes  was  15  m. 

The  source  shot  line  was  set  by  chaining  in  stakes  every  1.2  m  on  a  line  between  the 
two  boreholes.  Figure  30  is  a  plot  of  the  site  and  Fig.  23  shows  the  site  with  the  stakes 
installed.  Initial  source  spacing  was  derived  from  the  relations  given  in  Witten  and  Long 
(1986).  Using  a  representative  source  frequency  of  200  Hz  (from  Fig.  19)  and  assuming  an 
initial  sound  speed  in  the  soil  of  250  m/sec,  a  wavelength  of  1.25  m  was  calculated.  To  image 
an  object  with  a  radius  of  0.3  m  this  establishes  the  required  source  spacing  at  approximately 
0.6  m.  The  sources  were  set  to  provide  a  horizontal  coverage  of  from  1.2  to  13.4  m.  The 
maximum  source  distance  was  limited  by  available  working  space.  Depth  of  the  receiver 
array  (streamer)  was  selected  to  correspond  to  the  anticipated  depth  of  the  buried  pipe. 

The  field  station  for  these  tests  was  a  pickup  truck  outfitted  with  a  camper  top.  The  DAS 
and  supervisory  computer  were  operated  from  the  back  of  the  truck.  Power  was  supplied  by 
a  commercial  110- volt  650  watt  gasoline  powered  generator.  The  system  was  assembled  by 
wiring  the  computer  into  the  DAS,  the  hydrophones  to  the  individual  boards  in  the  DAS, 
and  the  accelerometer  on  the  BETSY  to  the  trigger  circuit  on  the  DAS.  The  streamer  was 
positioned  in  the  water  filled  borehole  by  setting  the  depth  of  the  first  hydrophone  to  a 
selected  depth  referenced  to  the  ground  surface. 

Since  this  test  was  the  first  designed  for  the  collection  of  data  to  support  GDT,  the  site 
was  selected  because  of  the  expected  simple  subsurface  conditions;  i.e.,  a  buried  pipe  in 
otherwise  homogeneous  clay  soil. 

Figure  31  is  a  gray-scale  image  of  relative  sound  speed  variations  derived  for  this  field 
experiment.  Most  important  is  the  clear  indication  of  the  pipe  centered  at  7.5  m  from  the 
receivers  and  1  m  deep.  Also  note  that  this  image  reveals  subsurface  conditions  more  complex 
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Figure  30:  Site  plan  for  the  Chestnut  Ridge  test  site. 
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CHESTNUT  RIDGE  SURVEY  -  MARCH  1987 


Figure  31:  Gray-scale  image  of  the  Chestnut  Ridge  survey  cross-section  from  the  March  1987 
field  study. 
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than  originally  believed.  The  elongated  shape  of  the  pipe  is  caused  by  the  offset  geometry  of 
the  source/receiver  configuration  as  discussed  in  Chapters  2  and  3  along  with  the  fact  that 
only  transmission  view  angles  were  considered.  Figure  32  is  a  synthetic  reconstruction  from 
perfect  computer  generated  data  of  a  circular  disk  in  a  homogeneous  media.  Comparison 
of  the  form  of  the  actual  result  with  the  form  of  the  optimum  achievable  result  for  the 
offset  vertical  seismic  profiling  configuration  clearly  demonstrates  that  GDT  was  successful 
in  finding  and  imaging  the  pipe.  The  other  features  identifiable  in  Fig.  31,  are  an  isolated 
area  of  high  sound  speed  located  close  to  the  receiver  array  (shown  as  a  black  area),  the 
shallow  horizontal  region  of  relatively  dry  soil  (shown  as  a  white  area),  and  a  transition 
with  depth  to  a  moist  soil  layer  (shown  as  light  gray)  resulting  from  an  extended  period  of 
rain  prior  to  the  field  survey.  Again,  the  system  configuration  tilts  these  inclusions  in  the 
direction  of  the  ray  path.  This  is  particularly  evident  for  the  layer. 

5.2  Chestnut  Ridge  October  1987  Survey 

Verification  of  the  methods  and  findings  from  the  initial  imaging  experiment  was  accom¬ 
plished  by  a  resurvey  of  the  site  in  October.  The  same  field  procedures  as  used  in  the  March 
survey  were  followed  here.  The  image  produced  from  this  study  is  presented  in  Fig.  33. 
This  study  was  conducted  during  extremely  dry  conditions  and  absent  are  the  fast  regions 
supporting  the  original  conclusion  of  their  presence  being  transient  elevated  soil  moisture. 
The  pipe  remains  in  approximately  the  same  position  as  seen  in  Fig.  31  as  does  the  fast 
shallow  area  near  the  borehole.  The  slight  differences  between  the  two  reconstructions  was 
caused  by  a  change  in  initial  receiver  positions.  For  the  first  study  the  top  receiver  was  0.6  m 
deep  while  in  the  second  it  was  only  0.15  m  deep. 

Also  evident  in  the  lower  left  of  this  image  is  an  isolated  region  of  elevated  sound  speed. 
This  is  likely  a  lens  of  more  compressed  soil.  Although  more  subtle  due  to  the  surrounding 
soil  moisture,  this  feature  does  appear  in  the  March  results  (Fig.  31). 

In  addition  to  verification  of  the  March  survey,  this  field  test  was  used  to  test  a  measure¬ 
ment  geometry  which  should  offer  improved  image  quality. 

One  image  artifact  that  can  occur  in  transmission  tomography  is  an  elongation  of  features 
in  the  predominant  direction  of  incident  wave  propagation.  This  is  a  result  of  the  limited 
propagation  or  projection  angles  inherent  in  most  measurement  geometries.  This  elongation 
results  from  the  resolution  always  being  superior  in  the  direction  normal  to  the  direction 
of  propagation  than  in  the  direction  parallel  to  the  direction  of  propagation.  Artifacts  of 
this  type  are  evident  in  some  of  the  previously  presented  images  and  can  be  clearly  seen  in 
Fig.  33  which  shows  an  image  of  a  cross-section  containing  a  0.61  m  diameter  buried  pipe 
and  other  identified  features.  This  image  is  derived  using  the  backpropagation  algorithm 
applied  to  data  collected  in  an  offset  VSP  configuration  with  a  borehole  containing  a  receiver 
array  located  to  the  left  of  the  imaged  cross-section.  For  this  geometry,  the  directions  of 
propagation  are  all  downward  from  right  to  left.  The  pipe,  which  should  display  a  circular 
cross-section,  is  clearly  elongated  in  the  direction  of  propagation.  A  means  to  minimize  this 
artifact  is  to  implement  a  measurement  geometry  which  offers  a  broader  range  of  detectable 
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Figure  32:  Gray-scale  image  of  a  circular  pipe  derived  from  simulated  data. 
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CHESTNUT  RIDGE  SURVEY  -  OCTOBER  1987 
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Figure  33:  Gray-scale  image  of  the  Chestnut  Ridge  survey  cross-section  from  the  October 
1987  field  study. 
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incident  propagation  angles.  Once  such  geometry  is  offset  VSP  using  receiver  arrays  in 
two  parallel  boreholes  with  sources  deployed  on  the  ground  surface  between  the  two  receiver 
arrays.  This  configuration  doubles  the  available  transmission  viewing  angles  that  can  be  used 
in  the  imaging  procedure.  This  configuration  was  implemented  for  the  same  cross-section 
as  shown  in  Fig.  33.  The  resulting  image  of  the  portion  of  the  cross-section  containing 
the  pipe  is  given  in  Fig.  34.  A  comparison  of  Figs.  33  and  34  demonstrates  a  considerable 
improvement  in  image  quality.  Slight  elongation  still  occurs;  however,  this  is  not  believed  to 
be  an  image  artifact  but  due  to  the  measurement  geometry  relative  to  the  pipe.  The  cross- 
section  defined  by  the  subsurface  region  between  the  two  boreholes  is  not  perpendicular  to 
the  axis  of  the  pipe.  In  this  geometry  the  cross-section  of  the  pipe  should  properly  have  a 
slight  elliptic  shape  as  portrayed  in  Fig.  34. 

5.3  Bear  Creek  Valley  March  1988  Survey 

Following  the  successful  completion  of  the  two  initial  field  studies,  it  was  appropriate  to 
conduct  a  more  controlled  field  study  directed  towards  the  imaging  of  buried  wastes.  This  was 
performed  in  the  Bear  Creek  Valley  area  west  of  the  Y-12  facility  on  the  ORNL  reservation. 
No  boreholes  were  developed  here,  but  rather,  an  existing  PVC-cased  monitoring  well  was 
used  to  house  the  hydrophone  array.  Coupling  was  maintained  by  periodically  adding  water 
to  the  weU.  Source  locations  were  established  in  the  same  manner  as  described  for  the 
Chestnut  Ridge  surveys. 

A  second  objective  of  this  field  study  was  to  determine  the  increased  data  acquisition 
speed  offered  by  the  COMPAQ  Portable  386  personal  computer  which  replaced  the  previously 
used  IBM-XT  just  prior  to  this  field  test.  Several  minor  software  changes  were  made  to  the 
data  acquisition  software  to  take  better  advantage  of  the  new  supervisory  computer. 

Two  types  of  waste-simulating  targets  were  selected  for  burial  at  predetermined  locations 
on  this  site.  These  targets  are  0.2  m^  metal  drums  approximately  0.61  m  in  diameter  and 
1  m  tall,  both  empty  and  water-filled,  as  well  as  plastic  bags  containing  styrofoam  packing 
peanuts  of  about  the  same  dimensions  as  the  metal  drums.  Six  targets  were  deployed  over 
the  two  subsurface  cross-sections  depicted  as  lines  A  and  B  in  Fig.  35.  An  array  of  29 
uniformly  spaced  hydrophones  spanning  a  0.61  to  4.9  m  depth  interval  was  located  in  a 
cased  monitoring  well.  Source  positions  were  established  every  0.61  m  over  two  survey  lines 
extending  radially  outward  from  the  well  containing  the  receiver  array  (Fig.  35).  Source 
positions  ranged  from  1.2  to  14.6  m  from  the  well  along  line  A  and  from  1.8  to  14.6  m  from 
the  well  along  line  B.  . 

The  resulting  images  of  these  two  cross-sections  are  displayed  in  Figs.  36  and  37.  The 
dominant  features  in  the  image  of  Line  A  are  a  fast  area  near  the  monitoring  well  identified  as 
an  area  of  folded  and  weathered  shale  observed  in  the  excavation  operations  having  a  sound 
speed  greater  than  the  surrounding  soil,  the  drum  and  the  packing  peanuts.  The  styrofoam 
packing  peanuts  appear  in  the  image  to  have  a  sound  speed  slightly  greater  than  the  local 
soil.  This  result  is  unexpected  because  of  the  high  compressibihty  of  the  target  and  may 
be  caused  by  a  blurring  between  this  target  and  the  fast  adjacent  weathered  shale.  Clearly 
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Figure  34:  Gray-scale  image  of  the  buried  pipe  cross-section  at  the  Chestnut  Ridge  site  using 
data  from  two  boreholes. 
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Figure  35:  Buried  target  configuration  at  the  two  cross-sections  developed  at  the  Bear  Creek 
site. 
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Figure  36:  Gray-scale  image  of  fast  (high  sound  speeds)  features  below  Line  A  at  the  Bear 
Creek  site. 
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Figure  37:  Gray-scale  image  of  cross-section  below  Line  B  at  the  Bear  Creek  site  for  (a)  full- 
range  plotting  contrast,  (b)  enhanced  and  inverted  contrast  to  emphasize  the  slow  (low  sound 
speed)  features,  and  (c)  enhanced  contrast  to  emphasize  the  fast  (high  sound  speed)  features. 


displayed  in  this  figure  as  black  at  about  1  m  deep  and  4.6  m  from  the  well  is  the  image  of 
the  two  water-filled  drums  which  exhibits  a  sound  speed  greater  than  both  the  surrounding 
soil  and  the  wedge  of  weathered  shale. 

Images  of  cross-section  B  are  given  in  Fig.  37.  Because  this  image  displays  a  number  of 
features,  all  of  which  cannot  be  clearly  shown  using  a  single  plotting  contrast,  results  are 
presented  for  three  plotting  contrasts.  The  first  [Fig.  37(a)]  is  a  full  range  plotting  contrast. 
All  buried  targets  appear  in  this  image  at  the  proper  location,  but  the  empty  drum,  as 
expected,  does  not  exhibit  a  strong  contrast  to  the  surrounding  soil.  The  low  sound  speed 
features  are  displayed  in  Fig.  37(b).  These  are  the  empty  drum  appearing  as  the  gray  area 
2  m  deep  and  9  m  from  the  well  and  the  bag  of  styrofoam  packing  peanuts  appearing  as 
the  white  area  1  m  deep  and  4  m  from  the  weU.  The  sound  speed  of  the  empty  drum  is 
slightly  greater  than  that  of  the  host  soil  while  the  sound  speed  of  the  styrofoam  peanuts  is 
slightly  less  than  that  of  the  soil.  The  image  of  the  water-filled  drum  is  less  blurred  when 
the  plotting  contrast  is  adjusted  to  highlight  the  high  sound  speed  features  [Fig.  37(c)]. 

5.4  Dinosaur  Site  Survey 

A  field  study  was  performed  in  March  1988  at  a  remote  site  in  the  high  desert  of  New  Mexico. 
The  objective  of  this  study  was  to  image  the  buried  skeletal  remains  of  a  large  sauropod 
dinosaur.  The  reasons  that  this  particular  project  was  undertaken  were  (1)  to  demonstrate 
that  the  system  could  operate  under  rugged  field  conditions  and  (2)  to  test  the  imaging 
method  in  a  different  geologic  setting,  sandstone  rather  than  the  clay  and  weathered  shale 
typical  of  the  ORNL  reservation.  While  all  other  field  tests  described  here  were  sponsored 
by  USATHMA,  this  task  was  completed  with  ORNL  discretionary  funds. 

No  boreholes  were  developed  at  this  site,  instead  existing  coreholes  were  used.  Capped 
PVC  pipes  were  placed  in  these  holes,  filled  with  water,  and  annular  regions  between  the 
hole  walls  and  the  casings  were  backfilled  with  sand  from  the  site.  The  BETSY  seismic  gun 
was  the  primary  energy  source  employed  at  the  site,  but  a  simple  hammer  and  steel  plate 
was  also  successfully  tested  here.  Figure  38  is  a  photograph  of  this  site. 

The  image  shown  in  Fig.  39  is  of  a  cross-section  of  the  skeletal  remains  of  the  sauropod 
dinosaur  buried  in  the  host  sandstone  formation.  The  black  area  in  this  figure  is  a  bone  from 
the  pelvic  region  of  the  dinosaur.  This  image,  which  was  developed  using  the  backpropa- 
gation  algorithm  in  an  offset  VSP  configuration,  is  included  to  demonstrate  an  important 
feature  of  remote  sensing  with  this  geometry.  Here,  the  subsurface  region  containing  the 
bone  mass  is  buried  beneath  an  area  of  numerous  surface  obstructions  (vegetation,  boulders, 
and  irregular  terrain).  At  this  location  it  Wcis  not  possible  to  deploy  sources  immediately 
above  the  desired  study  region.  Consequently,  source  locations  were  established  along  a  line 
radially  outward  from  the  borehole  containing  the  receiver  array  beginning  at  a  distance 
well  beyond  the  bone  location.  The  measurement  geometry  is  illustrated  in  Fig.  40.  The 
surface  features  at  this  site  would  make  detection  of  the  target  of  interest  quite  difficult  with 
any  remote  sensing  method  based  on  a  reflection  measurement  geometry;  however,  imag¬ 
ing  was  easily  accomplished  by  the  application  of  GDT.  This  result  is  significant  in  many 
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Figure  38;  Photograph  of  an  exposed  and  plastered  bone  mass  at  the  dinosaur  site. 
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Figure  39:  Gray-scale  image  of  a  vertical  cross-section  believed  to  contain  buried  dinosaur 
bone. 
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environmental  applications  in  that  it  demonstrates  the  “look-under”  capability  available  in 
tomography.  This  capability  is  important  not  only  for  naturally-occurring  surface  obstruc¬ 
tions  such  as  those  encountered  at  the  dinosaur  site,  but  for  man-made  obstructions,  such 
as  buildings,  above-ground  storage  tanks,  areas  of  surface  or  shallow  contamination,  etc.,  as 
well. 


5.5  Fort  Rucker  Survey 

The  most  recent  field  test  of  GDT  was  performed  at  Ft.  Rucker,  Alabama  during  December 
1988.  The  objective  of  the  effort  was  to  demonstrate  the  system  in  a  “production”  mode  at 
a  actual  waste  site.  To  achieve  this  goal  a  significant  upgrade  of  the  software  was  undertaken 
to  produce  the  operating  system  described  in  Chapter  4.  In  its  current  form,  the  system 
operates  in  a  mode  comparable  to  the  anticipated  commercial  technology. 

The  imaged  region  is  a  portion  of  a  trench  which  itself  is  part  of  a  series  of  trenches 
comprising  a  closed  landfill  at  Ft.  Rucker.  No  records  exist  of  the  contents  and  disposal 
practices  employed  during  its  operation,  the  landfill  was  closed  in  the  early  1970’s  and  is 
believed  to  predominantly  contain  wood  and  other  construction  debris.  In  order  to  identify 
areas  containing  isolated  inclusions  or  other  features  of  interest  in  a  background  setting  of 
largely  homogeneous  decaying  organic  materials,  a  site  screening  exercise  was  conducted 
approximately  one  month  prior  to  the  December  field  test.  At  this  time,  a  large  portion  of 
the  site  was  surveyed  using  electromagnetic  terrain  conductivity  mapping.  This  technique 
is  implemented  near  the  ground  surface  and  detects  subsurface  regions  of  elevated  electrical 
conductivity.  All  of  the  trench  areas  surveyed  exhibited  increased  conductivity  readings 
relative  to  background  due  to  numerous  small  metalhc  objects  distributed  within  each  trench. 
One  area  in  particular  exhibited  an  anomolously  high  metallic  content  and  this  location  was 
marked  and  the  area  surrounding  this  position  was  selected  as  the  primary  site  for  the  GDT 
field  test. 

Immediately  following  the  EM  survey,  the  site  was  prepared  for  the  application  of  GDT. 
This  consisted  of  minor  grading  mainly  to  provide  access  for  the  drill  rig  followed  by  the 
drilling  and  development  of  cased  boreholes.  Figure  41  is  a  photograph  of  the  site  after  site 
prep  work  was  completed. 

During  the  December  field  test  at  Ft.  Rucker,  data  was  collected  and  images  reconstructed 
for  five  vertical  cross-sections  utilizing  two  boreholes.  As  shown  in  Fig.  42,  three  lines  were 
surveyed  along  the  axis  of  a  trench  (Lines  2,  3,  and  4),  one  line  was  surveyed  across  the 
trench,  and  one  survey  line  was  established  away  from  the  trench  area  to  image  undisturbed 
background  conditions.  Pertinent  data  collection  parameters  for  each  of  these  cross-sections 
is  summarized  in  Table  1.  Figure  43  shows  the  horizontal  and  vertical  extents  of  the  three 
axial  trench  surveys.  Images  of  these  five  cross-sections  are  shown  in  Figs.  44-48.  It  should 
be  noted  that  in  all  these  images  a  full  range  of  view  angles  (transmission  and  backscatter) 
are  used  while  in  all  previous  results  only  transmission  angles  are  considered. 

Figure  44  is  an  image  of  Line  1  and  displays  the  background  conditions  typical  of  undis¬ 
turbed  areas  of  the  study  site.  Most  evident  in  this  image  are  numerous  clay  and  sand  soil 
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Figure  41:  Photograph  of  the  Ft.  Rucker  site. 
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Figure  42;  Site  plan  for  the  Ft.  Rucker  site. 
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Figure  43:  Illustration  of  the  three  imaged  cross-sections  along  the  axis  of  a  Ft.  Rucker 
landfill  trench. 
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Figure  44:  Gray-scale  image  of  the  cross-section  below  Line  1  (background  conditions)  at 
the  Ft.  Rucker  site. 
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Figure  45:  Gray-scale  image  of  a  cross-section  along  the  trench  axis  below  Line  2. 
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Figure  46:  Gray-scale  image  of  a  cross-section  along  the  trench  axis  below  Line  3. 
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Figure  47:  Gray-scale  image  of  a  cross-section  along  the  trench  axis  below  Line  4. 
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Figure  48:  Gray-scale  image  of  a  cross-section  below  Line  5  which  extends  across  the  trench. 
The  plotting  contrast  has  been  inverted  to  highlight  the  slow  (low  sound  speed)  features. 
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Table  1:  Data  acquisition  parameters  used  in  the  Ft.  Rucker  field  demonstration. 


Line  no. 

Minimum  receiver 
depth  (m) 

Minimum  source 
offset  (m) 

Source 
spacing  (m) 

Sampling 
interval  (//s) 

Frequency 

(Hz) 

1 

0.61 

1.21 

176 

2 

0.61 

1.21 

1.21 

98 

3 

2.44 

9.75 

176 

4 

0.30 

12.19 

1.21 

98 

5 

2.44 

1.21 

1.21 

100 

98 

strata  found  in  southeastern  Alabama.  Note  in  this  image  that  these  layers  exhibit  a  slight 
dip  which,  unlike  the  artifically  tilted  dry/moist  soil  contact  shown  in  Fig.  31  resulting  from 
geometric  constraints,  is  believed  to  be  real.  The  fastest  layer  occurs  at  the  surface  in  this 
image.  The  presence  of  this  layer,  extremely  dense  hardpan  soil,  and  its  dip  were  confirmed 
by  penetration  studies.  The  only  unexpected  feature  in  this  image  is  an  isolated  fast  inclu¬ 
sion  which  appears  as  the  black  area  in  the  upper  right  of  Fig.  44.  Excavation  at  this  precise 
position  revealed  a  small  mass  of  buried  concrete. 

Images  of  the  three  surveyed  cross-sections  along  the  trench  axis  are  provided  in  Figs.  45- 
47.  Figure  43  should  be  used  as  a  reference  when  comparing  these  results.  It  should  be  noted 
that  overlap  exists  among  these  three  cross-sections  and  that  in  these  overlapping  regions 
most  features  appear  consistently  within  individual  cross-sections.  It  is  also  important  to 
realize  that  each  plot  displays  a  full  range  of  gray  scales  determined  by  the  minimum  and 
maximum  sound  speed  occuring  within  that  particular  cross-section.  Therefore,  gray  levels 
are  image-specific  and  not  consistent  from  one  image  to  the  next.  For  this  reason  a  feature 
common  to  more  than  one  cross-section  may  appear  as  different  gray  scales. 

Below  Lines  2  and  4,  Figs.  45  and  47;  respectively,  the  mid-depths  of  the  trench  is 
characterized  by  very  slow  sound  speeds  displayed  in  these  two  images  as  white  or  light 
gray.  This  trench  is  believed  to  contain  wooden  construction  debris.  In  the  humid  climate 
of  Ft.  Rucker,  material  of  this  type  which  has  been  in  the  ground  for  almost  20  years, 
has,  no  doubt,  decayed  to  the  point  where  it  is  quite  compressible,  having  a  consistency 
similar  to  peat  moss.  Thus,  the  very  low  sound  speed  is  likely  to  be  a  result  of  the  high 
compressibility  of  buried  material  in  this  area.  In  both  of  these  figures,  the  highest  sound 
speed  areas  occur  near  the  top  and  bottom  of  the  images.  These  regions  represent  the  vertical 
boundaries  of  the  trench.  The  feature  appearing  as  black  in  the  lower  left  portion  of  the 
Line  4  image  (Fig.  47)  exhibits  the  greatest  sound  speed  of  all  the  Ft.  Rucker  trench  images. 
This  feature  is  approximately  4  m  long  and  its  position  corresponds  exactly  to  the  location 
of  the  anomolously  high  reading  observed  during  the  earlier  EM  survey.  This  would  suggest 
that  this  object,  or  objects,  is  metallic;  however,  the  reconstructed  sound  speed  is  not  nearly 
great  enough  for  this  to  be  a  solid  metal  mass.  This  again  is  consistent  with  the  EM  findings 
since  this  method  responds  to  metallic  surface  area  rather  than  metallic  mass.  This  feature 
could  be  a  large  number  of  drums  spaced  so  close  together  that  they  blur  together  in  the 
image,  one  or  several  larger  storage  tanks,  or  other  metallic  objects  with  large  surface  areas 
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but  low  metallic  mass. 

The  image  below  Line  3  (Fig.  46)  is  similar  to  the  two  other  trench  axis  figures,  but  is  too 
deep  to  display  the  upper  boundary  of  the  trench.  The  shallow  region  of  this  image  displays 
the  low  speed  typical  of  the  bulk  of  trench  contents  along  with  an  isolated  fast  inclusion 
in  the  upper  right.  Being  a  relatively  deep  cross-section,  most  of  the  area  displayed  in  the 
image  is  below  the  lower  boundary  of  the  trench.  The  only  feature  difficult  to  interpret 
in  the  image  is  the  slow  area  in  the  lower  left  corner  of  this  figure.  This  is  too  deep  to 
be  associated  with  the  landfill  operations.  It  could  be  a  naturally-occurring  feature,  but  is 
too  slow  for  this  to  be  a  reasonable  interpretation.  It  is  more  likely  that  this  is  an  image 
artifact  associated  with  either  edge  effects  of  the  slant  stack  or  the  fact  that  equipment 
failure  required  decommissioning  of  several  of  the  deep  receivers. 

The  final  image  from  the  Ft.  Rucker  field  study  is  the  survey  line  across  the  trench 
(Line  5).  This  cross-section  is  shown  in  Fig.  48.  Here,  the  plotting  contrast  has  been 
reversed  so  that  fast  areas  appear  as  lighter  gray  levels  and  slower  areas  appear  as  darker 
gray  levels.  In  this  figure,  the  trench  appears  as  the  dark  gray  and  black  area  in  the  upper 
central  portion  of  the  image.  As  in  the  trench  axis  figures  (Figs.  45-47)  the  trench  exhibits 
a  low  sound  speed  as  a  result  of  the  highly  compressible  fill  material.  This  image  clearly 
shows  the  bottom  and  lateral  trench  boundaries;  however,  the  imaged  region  spans  a  depth 
interval  too  deep  for  the  upper  trench  boundary  to  appear  in  this  figure. 
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6  CONCLUSIONS  AND  RECOMMENDATIONS 


This  report  presents  the  mathematical  basis  for  geophysical  diffraction  tomography  (GDT), 
a  description  of  the  specialized  field  instrumentation  and  its  implementation,  as  well  as  the 
results  of  all  field  studies  thus  far  performed.  These  results  suggests  that  GDT  can  fulfill  its 
intended  purpose,  i.e.  provide  high-resolution,  quantitative  images  of  subsurface  conditions 
at  existing  and  potential  hazardous  waste  storage  or  disposal  sites.  No  special  geologic 
criteria  were  specified  in  site  selection  to  improve  the  results.  A  known  waste  configuration 
represents  the  only  significant  fixed  variable  for  only  one  of  these  field  tests,  but  is  consistent 
with  the  overall  development  program  being  conducted.  The  field  tests  represent  an  iterative 
process  in  the  development  of  the  GDT  system.  In  the  first  effort  a  single  known  isolated 
inclusion  was  imaged.  Following  this,  experimental  conditions  have  progressed  to  analyses 
of  known  multiple  targets  of  both  greater  and  less  than  background  sound  speed  in  the  same 
imaging  plane,  and  finally  to  sites  with  no  a  priori  knowledge  of  subsurface  conditions. 

Recognizing  the  experimental  controls  applied,  the  results  obtained  provide  for  several 
very  important  conclusions  for  this  level  of  development.  Most  importantly,  all  of  the  buried 
objects  were  successfully  located  and  imaged.  Reconstructions  of  cross-sections  containing 
multiple  targets  characterized  by  both  weak  (sound  speeds  within  20%  of  background)  and 
strong  (sound  speeds  differing  from  background  by  an  order  of  magnitude)  inhomogeneities 
exhibited  image  quality  comparable  to  that  of  simple  cross-sections  containing  only  a  single 
target.  In  addition,  images  of  targets  having  a  sound  speed  less  than  background  were  of  com¬ 
parable  quality  to  those  having  a  sound  speed  greater  than  background.  Geologic  interfaces 
and  soil  moisture  were  not  a  problem  in  object  detection  and  identification,  but  represented 
additional  site  characterization  data  obtained  concurrently  in  the  imaging  process.  Finally, 
field  results  have  demonstrated  the  capability  to  detect  and  image  buried  objects  without 
sampling  needed  directly  over  the  target,  a  serious  constraint  with  existing  remote  sensing 
technology. 

In  the  most  recent  field  test  (Ft.  Rucker,  Alabama)  it  was  demonstrated  that  the  system 
could  be  operated  efficiently  in  a  “production”  mode.  In  this  test,  operations  were  conducted 
in  a  manner  similar  to  that  for  applications  at  actual  waste  sites.  Using  rapid  data  acquisition 
and  the  ability  to  immediately  display  results,  it  was  possible  to  design  the  field  program 
while  the  test  was  in  progress.  The  results  obtained  for  one  cross-sectional  survey  were  used 
in  establishing  subsequent  cross-sections  to  be  imaged.  In  this  production  mode,  data  could 
be  acquired  and  preprocessed  at  a  rate  better  than  one  source  firing  every  minute.  This  is 
comparable  to  the  time  required  to  reposition  the  source.  Because  of  the  level  of  automation 
programmed  into  the  system,  this  is  considerably  faster  than  conventional  seismic  methods. 
In  addition,  GDT  offers  much  higher  spatial  resolution  than  such  methods  along  with  the 
capability  to  display  final  images  in  the  field  rather  than  well  after  the  data  collection  has 
been  performed. 

Table  2  presents  the  time  required  to  survey  and  image  a  30  m  cross-section  as  a  function 
of  the  desired  smallest  resoluable  features.  This  table  was  developed  on  the  basis  of  the 
proven  resolution  of  one-quarter  of  a  wavelength,  a  required  half  wavelength  source  spacing. 
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and  a  demonstrated  capability  to  collect  data  at  a  rate  of  one  source  firing  per  minute  and  to 
perform  final  data  analysis  within  two  minutes.  It  can  be  seen  from  this  table  that  surveys 
can  be  performed  very  rapidly  and  that  survey  time  decreases  with  decreasing  resolution. 
Thus,  one  effective  way  to  quantify  areas  of  large  extent  is  to  first  operate  at  low  resolution 
to  establish  the  locations  of  disturbed  areeis,  such  as  trenches,  or  other  features  of  interest. 
Then,  based  on  these  results,  specific  sites  can  be  identified  for  higher  resolution  applications. 


Table  2:  Sampling  requirement  and  implementation  time  as  a  function  of  resolution  for  30  m 


survey  line. 


Resolution  (m) 

Source  spacing  (m) 

No.  of  sources 
(minutes) 

Survey  time 

0.3 

0.6 

50 

52 

0.6 

1.2 

25 

27 

1 

2 

15 

17 

1.5 

3 

10 

12 

3 

6 

5 

7 

5 

10 

3 

5 

While  the  existing  data  acquisition  system  has  proven  capable  of  efficient  field  opera¬ 
tions,  it  is  inadequate  as  a  prototype  for  a  commercially  viable  system.  This  is  because  it 
was  designed  for  establishing  the  data  acquisition  parameters  necessary  for  implementation 
of  GDT.  As  such,  it  was  designed  for  ea.se  of  modification  using  wire  wraps  rather  than  more 
rugged  printed  circuits.  This,  along  with  the  limited  computing  power  of  personal  computers 
and  microprocessor  chips  available  at  the  time  the  system  was  built,  made  it  unnecessarily 
large  and  cumbersome.  It  is  envisioned  that  the  next  generation  of  data  acquisition  system 
will  serve  as  a  prototype  for  a  commercial  product.  By  using  state-of-the-art  portable  per¬ 
sonal  computers  and  microprocessors,  the  large  data  acquisition  cabinet  would  be  replaced 
by  single  data  acquisition  board  installed  in  the  portable  supervisory  PC.  With  it,  the  en¬ 
tire  DAS  supervisory  computer  system  would  be  reduced  in  size  to  a  small  portable  PC. 
In  addition,  the  entire  system  would  be  battery-powered  eliminating  the  need  for  local  AC 
power  or  a  portable  generator.  The  awkward  streamer  containing  the  array  of  hydrophones 
and  preamps,  would  be  replaced  by  individual  hydrophone/preamp  pairs  which  could  be 
transported  to  the  field  in  a  container  the  size  of  a  large  briefcase.  The  receiver  array  would 
be  assembled  in  the  field  by  connecting  the  components  together  with  appropriate  lengths 
of  coaxial  cable. 

The  final  point  to  be  made  deals  with  demonstrated  and  anticipated  applications  of  GDT. 
The  method  is  capable  of  high  resolution  imaging  not  available  from  any  existing  remote 
sensing  technology.  Demonstrated  abilities  include: 

1.  detection  of  isolated  inclusions 

2.  location  of  interface  boundaries 
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3.  detection  of  elevated  moisture  content 

4.  identification  of  homogeneous  (clean)  areas. 

From  these  initial  results  numerous  direct  applications  can  be  inferred.  The  benefits  for 
buried  object  detection  are  obvious.  Spatial  delineation  of  waste  burial  areas  is  a  simple 
extension  of  imaging  isolated  inclusions.  Stratigraphy  and  identification  of  localized  water 
tables  are  two  uses  applicable  in  site  characterization  studies.  The  ability  to  locate  isolated 
soil  moisture  combined  with  the  look  under  capability  of  the  technique  offers  promise  in  iden¬ 
tification  of  leaking  ground-level  or  underground  storage  units.  While  all  images  presented 
in  the  previous  section  are  for  vertical  cross-sections,  Fig.  49  displays  a  sloping  cross-section 
which  represents  a  detection  surface  spanning  the  horizontal  extent  of  the  disposal  unit.  In 
this  manner  receivers  are  installed  at  a  constant  depth  on  one  side  of  the  storage  tank  and 
sources  are  fired  on  the  opposite  side  either  on  the  surface  or  downhole.  The  resultant  image 
could  detect  any  significant  areas  of  isolated  moist  soils  in  the  imaging  plane  which  would 
infer  leakage  from  the  tanks  or  ponds.  In  a  related  example,  an  imaging  plane  could  be  con¬ 
structed  parallel  to  a  subsurface  grout  curtain  or  other  vertical  barrier.  This  implementation 
would  allow  detection  of  all  significant  areas  of  soil  moisture  on  the  down  gradient  side  of 
the  barrier  thus  allowing  for  100%  evaluation  of  the  barrier  as  compared  to  the  less  certain 
method  of  well  sampling. 

Through  the  results  and  potential  applications  presented  here,  it  is  seen  that  innovative 
implementation  of  the  proven  capabilities  of  geophysical  tomography  does  suggest  numerous 
opportunities  for  use  in  environmental  engineering.  Further  development  and  applications 
will  serve  to  further  define  and  improve  the  capabilities  of  the  technology;  however,  the 
system  needed  is  ready  for  its  primary  use,  characterization  of  hazardous  waste  burial  sites. 
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Figure  49;  Illustration  of  a  geophysical  diffraction  tomography  measurement  configuration 
which  can  be  used  to  look  under  liquid  waste  storage  or  disposal  units.  Such  a  geometry 
could  be  used  for  the  early  identification  of  leakage. 
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A  MATHEMATICAL  CONSIDERATIONS 

A.l  Inversion  Procedure 

The  propagation  of  acoustic  waves  is  governed  by  the  scalar  wave  equation 

where  p{x_,t)  is  the  pressure  as  a  function  of  position  x  and  time  t,  c{x)  is  the  spatially- 
variable  wave  speed,  and  f{x,t)  characterizes  the  spatial  and  temporal  distribution  of  the 
energy  source.  It  is  more  convenient  to  consider  the  frequency  dependence  than  the  time 
dependence,  so  that  eq.  (A.l)  is  Fourier  transformed  in  time  to  yield 


V^p  +  kQn^{x)p  =  f{x,Lo) 

(A.2) 

where 

p(x,t<;)  =  j  p{x,t) 

(A.3a) 

=  J  J{x,t) 

(A.3b) 

n(x)  =  Co/c(x)  is  the  spatially-variable  refractive  index,  Cq  is  a  reference  sound  speed, 
u  is  the  angular  frequency  (radians/sec),  and  kQ  =  is  the  reference  wavenumber  at 

frequency  lo. 

The  imaging  problem  is  now  one  of  extracting  information  about  n(x)  from  measured 
values  of  p.  The  most  rigorous  way  to  accomplish  this  is  by  direct  inversion  of  eq.  (A. 2). 
This  approach  is  not  practical,  so  it  becomes  necessary  to  linearize  this  equation  by  invoking 
a  weak  scatter  approximation.  This  is  done  by  first  expressing  the  refractive  index  as 

”^(2.)  =  1  —  e0(2.)  (A. 4) 

where  0{x^  is  referred  to  as  the  “object  profile”  and  e  is  the  small  parameter  to  explicitly 
represent  a  weak  scatter  approximation.  For  homogeneous  conditions,  n{x)  =  1  and 

0{x)  =  0.  Thus,  for  weak  scatter  situations,  the  object  profile  is  small. 

The  first  type  of  weak  scatter  approximation  to  be  discussed  is  the  Born  approximation. 
Within  this  approximation,  it  is  assumed  that  the  pressure  p  is  composed  of  the  sum  of 
the  pressure  resulting  from  propagation  in  homogeneous  conditions,  this  is  referred  to  as  the 
incident  pressure  pi,  and  a  perturbed  pressure,  p*,  associated  with  nonzero  values  of  0{x) 
(inhomogeneities).  This  is  written  as 

p  =  p.  +  eps  (A. 5) 

where  the  small  parameter  e  is  explicitly  included  to  reflect  the  fact  that  only  a  small 
perturbation  in  pressure  will  result  from  inhomogeneities  of  the  order  eO{£).  Substitution 
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of  eqs.  (A. 4)  and  (A. 5)  into  eq.  (A. 2),  equating  terms  of  equal  powers  of  e  and  neglecting 
terms  of  order  yields  the  set  of  linearized  equations 


V^p,  +  kip,  =  f 


(A.6) 


and 

V^Ps  +  kip,  =  kl  0{x)pi,  (A.7) 

where  eq.  (A.6)  governs  propagation  in  homogeneous  conditions  [0(x)  =  0],  while  eq.  (A.7) 
governs  the  corrective  term  for  the  pressure  due  to  inhomogeneities. 

Another  weak  scatter  approximation  is  the  Rytov  approximation.  Here,  the  total  pressure 
p  is  represented  by 

p  ^  ^  p.e^^'  =  ,  (A.8) 


where  tp  is  referred  to  as  the  complex  phase,  ipo  is  the  incident  or  unperturbed  phase 
(p-  =  e’^®)  and  ip'  is  the  perturbed  phase  associated  with  nonzero  values  of  0{x).  The  real 
part  of  Ip  contains  amplitude  information  while  the  imaginary  part  represents  variations  in 
wave  phase.  Substituting  eq.  (A.8)  into  eq.  (A. 2)  and  collecting  terms  of  equal  powers  of  e 
yields 


(VVi)e'’^'  +  2e(Vp,)  •  {ViP')  +  ep,-  ( W) 

+  kl  [  1  -  eO{x)  ]  Pi  =  fU,u;)  .  (A.9) 

Utilizing  the  Rytov  transform  D  —  p,t/’',  eq.  (A.6)  and  the  relationship  =  I  —  eip' 

for  eip'  small,  eq.  (A.9)  simplifies  to 

+  klD  =  klO{x)p,  .  (A.IO) 


The  two  weak  scatter  approximations  can  be  unified  by  defining  D  to  be  the  reduced 
data,  where 


D{x,u:) 


Piip'  ,  p  =  pic'^'^  ;  Rytov  approximation 
Ps  ,  P  =  Pi  +  Ps  Born  approximation 


(A.ll) 


and  both  satisfy  the  linear  equations 


V^D  +  klD  =  kl  0{x)p,  ,  (A.12) 

VVi  +  klpi  =  .  (A. 13) 

Following  is  the  derivation  of  an  analytic  inversion  procedure  for  the  above  system  of 
equations.  This  is  the  mathematical  basis  for  diffraction  tomography.  The  measurement 
geometry  considered  here  is  the  relatively  general  one  shown  in  Fig.  A-1.  Both  sources  and 
receivers  are  deployed  along  lines  having  arbitrary  orientation  both  to  each  other  and  the 
horizontal.  The  sources  are  located  at  positions  £o  on  the  line  defined  by  unit  vectors  n 
and  212  while  receivers  are  positioned  at  £'  along  a  line  with  unit  vectors  rn  and  22l2- 
The  absolute  position  relative  to  the  reference  coordinates  x,y  of  a  source  and  receiver  are 


84 


ORNL-DWG  87-18524 


Figure  A-1:  Geometry  and  notation  used  in  the  derivation  of  the  inversion  algorithm. 
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defined  by  the  vectors  and  r',  respectively.  The  solution  to  eq.  (A.  12)  for  measured 
values  of  reduced  data  D  at  position  r'  is  given  by  (Morse  and  Feshbach  1953) 


D(l')  =  -'y  j<‘L  Pi()  0({)  Ho  (*l£'  -  £ I)  (A.14) 

where  Hq  is  the  zeroth  order  Hankel  function  of  the  first  kind.  This  function  can  be 
expressed  in  terms  of  its  angular  spectrum  by 


Hoiko  |x|) 


da 


e 


i(  am2 


v/*o  -  Ja) 


X 


(A.15) 


Substituting  eq.  (A.15)  into  eq.  (A.14)  yields 

+  y/*'o  ~  ai)  z' 


■  j  di  pii)  0(£)e-’<“"‘^  +  la)  •  i  .  (A.16) 

The  problem  of  imaging  is  now  one  of  inverting  eq.  (A.16)  to  express  the  desired  quantity, 
0(x),  in  terms  of  the  known  values  of  D.  This  is  most  easily  accomplished  for  plane  wave 
illumination  and,  although  plane  waves  cannot  be  generated  in  the  field;  for  illustration 
purposes,  this  procedure  is  outlined  below. 

For  a  plane  wave  having  a  propagation  direction  defined  by  the  unit  vector  the 

incident  field  is  simply 

p,{x)  =  ,  (A.17) 

so  that  D  for  plane  wave  illumination,  denoted  as  Dp^,  can  be  written  as 


da 


\A^- 


0“ 


^  /-T^ 

47r  J  _  ^2 

•  j  d^  (9(^)e-'(om2  +  yA|TrT^  m  -  fcolo)  •  i 

By  defining  the  spatial  Fourier  transform  operator 

Dpw  =  I  dt'  ■  t'  Dp^{r:) 


i(am2  +  _  (,2 


(A.18) 


(A.19) 


and  applying  this  operator  to  both  sides  of  eq.  (A.19)  yields 


7p  ~  Sk  ■  l'  ^ 

DpW  (r'  ■  rn,K)  = - - - .  -0[  ko{s  -  ^)  ]  ,  (A. 20) 

^  V^o  - 
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where  5  =  ^  |  ^/k'^  —  k?  m  +  «:m2 

Equation  (A. 20)  is  the  holographic  lens  equation  and  relates  the  one-dimensional  spatial 
Fourier  of  the  data  along  the  receiver  line  to  the  two-dimensional  spatial  Fourier  transform 

0  [  koU  -  ^)]  =  jdi  Oii)  ^  (A.21) 


of  the  object  profile. 

In  real  applications,  acoustic  sources  are  typically  impulsive  point  sources.  The  incident 
field,  Pi,  at  a  position  x  for  a  point  source  located  at  point  in  two-dimensions  is 
governed  by 

Pi  +  koPt  =  -  A-k  6{x  -  r^)  (A.22) 

and  the  solution  to  this  equation  is 


p.(^)  =  iTT  Ho  {  ko  \  X  -  \)  = 

^«(on2  -(-  yjkl  -  q'^  2.)  •  -  ^o) 

Ao  -  f 


(A.23) 


Substituting  this  expression  for  the  incident  field  into  eq.  (A.  16)  yields  the  expression  for 
the  measured  data  D  for  an  impulsive  point  source 


D{r')  =  —  f  f  _ dadq _  ^  ^i{Qm2  +  y/k^-a^m)-r'-i{qn2  +  y/kl-^n)  -r^  24) 

J  ^ki  -  -  ,2 

■  J  d^  0(0 


This  expression  can  be  simplified  by  defining  the  operator 

D{r',  i')  =  j  dio  0(2:')  • 


(A.25) 


Applying  this  operator  to  both  sides  of  eq.  (A. 26)  and  using  the  relationship 


j  dl  e‘^  ~  “  9^  2  )  •  ro 

=  27r  8{v  -  9)e-V^^r^  s  •  lo  (A.26) 


yields 


e-V^o  ^-^-o  r 

„2  J 


da 


J{  ocm2  +  \/kl  -  of2  m  )  •  r' 


\/ko  -  1/2  ^  \]kl  -  a2 

J  d^  0(0  e”*^  \/*'o  “  21  -  ^n.2  -\/*^o  -  2  )  • 


(A.27) 
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By  defining  ko  ■  ^  =  i/n2  +  ^k^  —  n, 


?  -  /  2^  V  u 

•  j  Moe"-^ 

and  applying  the  operator  defined  by  eq.  (A. 19),  eq.  (A. 27)  becomes 

iko  e-«^o  -  m  •  ii' 


g-'V^O  -  s.  ■  Lo 

(A.28) 

D  , 

(A.29) 

Dp^  (r'  •  m,  k)  =  - 


\jkl  - 


0\  ko  (s  -  ^)  ] 


(A.30) 


This  is  the  equation  of  a  holographic  lens  identical  to  that  for  plane  wave  illumination 
[  eq.  (A. 20)  ].  Thus,  the  operator  defined  by  eq.  (A. 20)  is  a  synthetic  aperture  lens  which 
produces  the  response  for  plane  wave  illumination  at  an  insonifying  angle  ^  given  by 
eq.  (A. 28). 

The  final  step  of  the  inversion  procedure  is  to  recover  0{x)  from  eq.  (A.30).  This  is 


- 


ni  ■  T 


■  b  e 


(A,31) 


where 

K.  —  ko{s  —  Sx))  .  (A. 32) 

For  an  offset  VSP  configuration,  n  =  (0,-1),  22,2  =  (1,0),  m  =  (1,0),  m2  =  (0,1)  and  by 
defining  ^  =  (cos0,  sin0),  the  wavevector  K_  can  be  expressed  in  terms  of  /c  and  6  by 


K  =  {K,,  Ky)  = 


\Ai 


— 


—  kocos9,  K  —  ko  sin0 


(A.33) 


Computationally,  it  is  more  practical  to  integrate  over  6  and  k  than  over  Kx  and  Ky. 
To  accomplish  this,  the  change  of  variables 


dKxdKy 


j 

dh\ 

dK-r 

dn 

de 

dK^ 

dKy 

dKdO  = 

df\ 

as 

—kosinO 
ko  cos  0 


dndO 


dndO  (A. 34) 

is  introduced.  For  the  offset  VSP  configuration,  realizable  viewing  angles  are  0  <  0  <  —  tt. 
For  this  range  of  viewing  angles,  the  mapping  defined  by  eq.  (A. 34)  is  not  injective  when 
(1)  — 7r/2  <  ^  <  0  and  k  —  kosmO,  and  (2)  — tt  <  6  <  — 7r/2  and  k  =  —  kosinO.  In 


sin  0  —  K  cos  0 
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addition,  it  is  necessary  to  remove  double  coverage  which  occurs  when  K^,  Ky  <  0  as  well 
as  evanescent  modes  |k|  >  ko.  This  is  accomplished  by  defining  the  filter  function 


F{K,e) 


0,  for  |k1  >  ko; 

0,  for  —  ko  <  K  <  ko  sin0  and  —  |  ^  <  0; 

0,  for  ~  k  s'm6  <  K  <  ko  and  —  tt  <  0  |  ; 

2  \\JkQ  —  K?  sin0  —  Kcos^l  ?  otherwise. 


with  this  the  final  form  of  the  expression  for  0(x)  becomes 


0(x) 


Iff  -ikolx cose  +  ysinS) 

2irHo  ' 

^2  g’V^^o  “  ~  51-ri)  +  i>iy 


(A.35) 


(A.36) 


A. 2  Phase  Unwrapping 

Implementation  of  geohysical  diffraction  tomography  within  the  Rytov  approximation  re¬ 
quires  a  quantification  of  the  perturbed  phase,  t/j'  [eq.  (A.  11)],  the  imaginary  part  of  which 
is  essentially  the  difference  in  distance,  measured  in  radians  or  wave  cycles,  between  the  full 
wave  field  and  the  wave  field  which  would  propagate  in  homogeneous  conditions.  For  an 
impulsive  source,  the  Fourier  transformed  pressure  will  be  complex 

p  ^  Pr  +  tPi  (A.37) 


having  real  and  imaginary  parts  Pr  and  Pi,  respectively.  Phase  perturbations  must  be 
deduced  from  this  relation  and  eq.  (A. 8)  where 


Pfi  +  iPj  = 

(A.38) 

The  real  part  of  ip,  ipn,  where  xp  =  '1^11  +  ixpj  is  simply 

=  tn{PPhTTj)  , 

(A.39) 

while  the  principal  value  of  the  imaginary  part  of  %p,  ip]f\  is 

given  by 

=  arctan  [Pij Pr)  . 

(A.40) 

The  arctan  only  returns  values  between  —it  and  tt,  the  principal  value,  which  represents 
only  the  fraction  portion  of  the  total  number  of  wave  cycles.  Thus,  xpi  is  related  to  tp\’’^ 

by 

+  27rn  (A. 41) 

where  n  =  0,  ±1,  ±  2,... 

The  phase  unwrapping  procedure  is  a  logical  method  for  deducing  n  based  on  available 
information.  To  accomplish  this,  two  assumptions  are  made.  First,  it  is  assumed  that  n 


89 


is  known,  based  on  field  information  for  the  shallowest  receiver  and  second  that  changes  in 
phase  are  smooth  with  receiver  depth  so  that  any  large  change  in  which  occurs  suggests 
a  point  of  phase  unwrap.  By  testing  the  difference  between  successive  values  in  depth  of 
points  of  large  change  in  the  difference,  say  ±  7r/2,  can  be  identified.  At  all  such 
points  a  value  of  it  27r  is  added  to  along  with  a  value  of  ±  2nn  associated  with  the 
xpj  from  the  previous  receiver  position. 

Having  accomplished  the  phase  unwrap  by  this  procedure,  the  phase  perturbation  is 
easily  computed  by  subtracting  the  complex  quantity  where 

=  p.  =  *7r  Ho  {koR)  ,  (A.42) 


from  the  full  complex  phase  xj;. 

A. 3  Relationship  of  Sound  Speed  to  Mechanical  Properties 

For  an  inviscid  fluid,  the  equation  governing  the  perturbed  pressure  D  [eq.  (A. 11)]  is 

V^£>  +  klD  =  -ki  I  ^,{x)  +  7.(1)  )  p,  -  1  (7.)  PI  ,  (A.43) 

where 

IpU)  =  [  /?(r)  -  00  ]/0o  (A.44) 

is  the  normalized  compressibility  perturbations  and 

lp{x)  =  [  pix)  -  Po  ]/p{x)  (A.45) 

is  the  normalized  density  perturbations.  Comparing  eq.  (A.43)  to  eq.  (A. 10)  it  is  easily  seen 
that 

0{x)  =  1  -  cl/c^{x)  =  -  [  7^(x)  +  7p(x)  \  7pU)  •  (A.46) 

This  clearly  shows  that  wave  speed  variations  are  a  result  of  perturbations  in  the  more 
fundamental  mechanical  properties,  density  and  compressibility.  A  similar  relationship  can 
be  obtained  from  the  acoustic  approximation  of  the  elastic  wave  equation.  Here  sound  speed 
is  a  composite  function  of  density  and  bulk  modulus. 
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